##############
# Montecarlo #
##############
library("MASS")
library("expm")
library("parallel")
#library("expm"")

cl <- makeCluster(detectCores() - 2)


# FUNCTIONS
{
eigenvector_assumption = function(Psi) # Qui si inserisce l'ussumption su autovettori
{
  r = ncol(Psi)
  # out = rep(1, r)       # lascio Psi com'e'
  # out = sign(Psi[1,])   # Primo elemento di ogni colonna positivo
   out = sign(apply(Psi, 2, mean)) # media di ogni eigenvector positiva
  
  return(out)
}
population_quantities = function(factors0)
{
  TT = nrow(factors0)
  r = ncol(factors0)
  
  sigma2 = r
  Sigma_lambda = diag(1, r)
  
  Sigma_F_hat = t(factors0) %*% factors0 / TT
  supp = sqrtm(Sigma_lambda) %*% Sigma_F_hat %*% sqrtm(Sigma_lambda)
  supp = eigen(supp)
  U = diag(supp$values, r)
  U_inv = solve(U)
  Psi = supp$vectors
  L = eigenvector_assumption(Psi)           
  Psi = Psi %*% diag(L)
  
  Q = sqrtm(U) %*% t(Psi) %*% solve(sqrtm(Sigma_lambda))
  H = sqrtm(Sigma_lambda) %*% Psi %*% sqrtm(U_inv) 
  
  supp = 3 + 27/(TT^3)*sum(diag((factors0 %*% t(factors0)))^2) + 18*r/TT
  sigma4 = as.numeric(solve(supp) * 3)
  
  out = list(U_inv, Psi, Q, H, sigma2, sigma4, Sigma_lambda, Sigma_F_hat)
  return(out)
}

# Factors
OLS = function(X, k)
{
  TT = nrow(X)
  N = ncol(X)
  factor = sqrt(TT) * as.matrix(eigen(X %*% t(X))$vectors[,1:k])
  loadings = 1/TT *(t(X) %*% factor)
  out = 1/(N*TT) * sum((X - factor %*% t(loadings))^2)
  return(out)
}
nfac_selectioncriteria = function(X, k, eps)
{
  eps1 = eps[1]
  eps2 = eps[2]
  TT = nrow(X)
  N = ncol(X)
  out1 = OLS(X, k)
  out3 = k * ((log(N))^eps1 * N^eps2)/sqrt(N)
  NEW1 = out1 + out3
  return(c(NEW1, out1))
}
nfac_optimization = function(data, kmax, eps, std = F, est_extract = TRUE)
{
  TT = nrow(data)
  X_temp = as.matrix(data)
  if (std == TRUE)
  {
    m = apply(X_temp, 2, mean)
    s = apply(X_temp, 2, sd)
    X_temp = (X_temp - rep(1, TT) %*% t(as.matrix(m)))/(rep(1, TT) %*% t(as.matrix(s)))
  }
  
  if (est_extract == T)
  {
    out = extraction(X_temp, kmax)
    return(out)
  }
  
  temp = matrix(NA, nrow = kmax, ncol = 2)
  for (k in 1:kmax) temp[k,] = nfac_selectioncriteria(X_temp, k, eps)
  out = rep(NA, 2)
  for (l in 1:2) out[l] = which(temp[,l] == min(temp[,l]))
  return(out)
}
extraction = function(X, nk)
{
  TT = nrow(X)
  N = ncol(X)
  
  eig = eigen(1/(TT*N) * X %*% t(X))$values[1:nk]
  
  val = rep(NA, nk)
  m_loadings = matrix(0, ncol = nk, nrow = nk)
  s_loadings = matrix(NA, ncol = nk, nrow = nk)
  factors = matrix(NA, ncol = nk, nrow = TT)
  sigma_4 = rep(NA, nk)
  gamma_twopass = matrix(NA, ncol = nk, nrow = nk)
  
  for (k in 1:nk)
  {
    temp = sqrt(TT) * as.matrix(eigen(1/(TT*N) * X %*% t(X))$vectors[,1:k])
    loadings = 1/TT *(t(X) %*% temp)
    m_loadings[1:k, k] = apply(loadings, 2, mean)
    val[k] = 1/(N*TT) * sum((X - temp %*% t(loadings))^2)
    temp1 = sum(diag((temp %*% t(temp))/TT)^2)
    temp1 = 3 + 27/TT*temp1 + 18*k/TT
    temp2 = 1/(N*TT) * sum((X - temp %*% t(loadings))^4)
    sigma_4[k] = temp2/temp1
    gamma_twopass[1:k,k] = solve(t(loadings) %*% loadings) %*% t(loadings) %*% apply(X, 2, mean)
  }
  factors[1:TT, 1:k] = temp
  s_loadings[1:k, 1:k] = (t(loadings) %*% loadings)/N
  
  out = list(val, eig, m_loadings, s_loadings, factors, sigma_4, gamma_twopass)
  names(out) = c("value_function", "eigenvalues", "mean_loadings", "sigma_loadings", "factors", "sigma_4", "gamma_twopass")
  return(out)
}

Ftilde_acm0_fun = function(t, factors0) # true par
{
  r = ncol(factors0)
  TT = nrow(factors0)
  
  Imatr = diag(1, r)

  temp = population_quantities(factors0)
  U_inv = temp[[1]]
  Psi = temp[[2]]
  Q = temp[[3]]
  H = temp[[4]]
  sigma2 = temp[[5]]
  sigma4 = temp[[6]]
  Sigma_lambda = temp[[7]]
  
  factors0 = matrix(factors0[t,], ncol = 1)
  
  A = rbind(Imatr, Imatr, Imatr)
  
  B11 = sigma4/TT * U_inv^2 + sigma4/(TT^2) * U_inv %*% t(H) %*% factors0 %*% t(factors0) %*% H %*% U_inv
  B12 = diag(0, r)
  B13 = diag(0, r)
  B21 = diag(0, r)
  B22 = sigma2 * U_inv
  B23 = sigma2/TT * U_inv %*% Q  %*% Sigma_lambda %*% factors0 %*% t(factors0) %*% H %*% U_inv
  B31 = diag(0, r)
  B32 = t(B23)
  B33 = c(sigma2/TT * (t(factors0) %*% Sigma_lambda %*% factors0)) * U_inv^2
  B1 = cbind(B11, B12, B13)
  B2 = cbind(B21, B22, B23)
  B3 = cbind(B31, B32, B33)
  B = rbind(B1, B2, B3)
  
  out = t(A) %*% B %*% A
  return(out)
}
Ftilde_acm_fun = function(t, k, factors, eigenvalue, valuef, sigma4) # est 
{
  TT = nrow(factors)

  Imatk = diag(1, k)
  sigma2 = TT/(TT - k)*valuef
  Sigma_lambda = diag(eigenvalue, k) - sigma2/TT*Imatk
  U_inv = solve(Sigma_lambda)
  
  factors = matrix(factors[t,], ncol = 1)
  
  A = rbind(Imatk, Imatk, Imatk)
  B11 = sigma4/TT * (U_inv %*% (Imatk + 1/TT * (factors %*% t(factors))) %*%  U_inv)
  B22 = sigma2 * U_inv
  B23 = sigma2/TT * U_inv %*% Sigma_lambda %*% factors %*% t(factors) %*% U_inv
  B32 = t(B23)
  B33 = c(sigma2/TT * (t(factors) %*% Sigma_lambda %*% factors)) * U_inv^2
  B12 = B13 = B21 = B31 =  diag(0, k)
  B1 = cbind(B11, B12, B13)
  B2 = cbind(B21, B22, B23)
  B3 = cbind(B31, B32, B33)
  B = rbind(B1, B2, B3)
  
  out = t(A) %*% B %*% A
  return(out)
}

Ftilde_acm_funbai = function(t, k, factors, eigenvalue, valuef, sigma4) # est 
{
  TT = nrow(factors)
  
  Imatk = diag(1, k)
  sigma2 = TT/(TT - k)*valuef
  Sigma_lambda = diag(eigenvalue, k) - sigma2/TT*Imatk
  U_inv = solve(Sigma_lambda)
  
  factors = matrix(factors[t,], ncol = 1)
  
  A = rbind(Imatk, Imatk, Imatk)
  B11 = sigma4/TT * (U_inv %*% (Imatk + 1/TT * (factors %*% t(factors))) %*%  U_inv)
  B22 = sigma2 * U_inv
  B23 = sigma2/TT * U_inv %*% Sigma_lambda %*% factors %*% t(factors) %*% U_inv
  B32 = t(B23)
  B33 = c(sigma2/TT * (t(factors) %*% Sigma_lambda %*% factors)) * U_inv^2
  B12 = B13 = B21 = B31 =  diag(0, k)
  B1 = cbind(B11, B12, B13)
  B2 = cbind(B21, B22, B23)
  B3 = cbind(B31, B32, B33)
  B = rbind(B1, B2, B3)
  
 # out = B22 #t(A) %*% B %*% A
  return(B22)
}



# Risk Premia
RiskPremiaPost_fun = function(factors0)
{
  out = apply(factors0, 2, mean)
  return(out)
}
RiskPremiatilde_fun = function(k, factors)
{
  factors = as.matrix(factors[,1:k])
  out = apply(factors, 2, mean)
  return(out)
}

RiskPremiatilde_acm0_fun = function(factors0)
{
  TT = nrow(factors0)
  r = ncol(factors0)
  
  Imatk = diag(1, r)
  
  gamma = RiskPremiaPost_fun(factors0)

  temp = population_quantities(factors0)
  U_inv = temp[[1]]
  Psi = temp[[2]]
  Q = temp[[3]]
  H = temp[[4]]
  sigma2 = temp[[5]]
  sigma4 = temp[[6]]
  Sigma_lambda = temp[[7]]
  
  A = rbind(Imatk, Imatk, Imatk)
  C11 = sigma4/(TT^2)*U_inv^2 + sigma4/(TT^2)*U_inv %*% t(H) %*% gamma %*% t(gamma) %*% H %*% U_inv
  C22 = sigma2/TT*U_inv
  C23 = sigma2/TT*U_inv %*% Q %*% Sigma_lambda %*% gamma %*% t(gamma) %*% H %*% U_inv
  C32 = t(C23)
  C33 = sigma2/TT*c((t(gamma) %*% Sigma_lambda %*% gamma))*U_inv^2
  C12 = C13 = C21 = C31 = diag(0, r)
  C1 = cbind(C11, C12, C13)
  C2 = cbind(C21, C22, C23)
  C3 = cbind(C31, C32, C33)
  C = rbind(C1, C2, C3)
  
  out = t(A) %*% C %*% A
  return(out)
}
RiskPremiatilde_acm_fun = function(k, factors, eigenvalue, valuef, sigma4)
{
  TT = nrow(factors)

  Imatk = diag(1, k)
  sigma2 = TT/(TT - k)*valuef
  Sigma_lambda = diag(eigenvalue, k) - sigma2/TT*Imatk
  U_inv = solve(Sigma_lambda)
  
  gamma = RiskPremiatilde_fun(k, factors)
  
  A = rbind(Imatk, Imatk, Imatk)
  C11 = sigma4/TT^2*U_inv^2 + sigma4/TT^2*U_inv %*% gamma %*% t(gamma) %*% U_inv
  C22 = sigma2/TT*U_inv
  C23 = sigma2/TT*U_inv %*% Sigma_lambda %*% gamma %*% t(gamma) %*% U_inv
  C32 = t(C23)
  C33 = sigma2/TT*c((t(gamma) %*% Sigma_lambda %*% gamma))*U_inv^2
  C12 = C13 = C21 = C31 = diag(0, k)
  C1 = cbind(C11, C12, C13)
  C2 = cbind(C21, C22, C23)
  C3 = cbind(C31, C32, C33)
  C = rbind(C1, C2, C3)
  
  out = t(A) %*% C %*% A
  return(out)
}

# SDF
SDFPost_fun = function(factors0)
{
  TT = nrow(factors0)
  r = ncol(factors0)
  
  rf = rep(1, TT)
  
  f_bar = apply(factors0, 2, mean)
  
  IvecT = matrix(1, nrow = TT, ncol = 1)
  ImatT = diag(1, TT)
  M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
  omega_hat = 1/TT * (t(factors0) %*% M_1T %*% factors0)
  omega_hat_inv = solve(omega_hat)
  
  out = 1/rf - 1/rf*t(f_bar) %*% omega_hat_inv %*% t(factors0 - rep(1, TT) %*% t(f_bar))
  return(c(out))
}
SDFtilde_fun = function(k, factors)
{
  TT = nrow(factors)

  rf = rep(1, TT)
  
  factors = as.matrix(factors[,1:k])
  
  gamma = apply(factors, 2, mean)
  
  omega = diag(1, k) - gamma %*% t(gamma)
  omega_inv = solve(omega)
  out =  1/rf - 1/rf*t(gamma) %*% omega_inv %*% t(factors - rep(1, TT) %*% t(gamma))
  
  return(c(out))
}

commutation_matrix = function(r, c)
{
  er = diag(1, r)
  ec = diag(1, c)
  
  out = matrix(0, ncol = r*c, nrow = r*c)
  for (i in 1:r)
    for (j in 1:c)
    {
      e1 = as.matrix(er[,i])
      e2 = as.matrix(ec[,j])
      out = out + (e1 %*% t(e2)) %x% (e2 %*% t(e1))
    }
  return(out)
}
SDFtilde_acm0_fun = function(t, factors0)
{
  TT = nrow(factors0)
  r = ncol(factors0)
  
  rf = rep(1, TT)
  
  IvecT = matrix(1, ncol = 1, nrow = TT)/TT
  ImatT = diag(1, TT)
  Imatk = diag(1, r)
  M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
  K_rT = commutation_matrix(r, TT)
  K_Tr = commutation_matrix(TT, r)
  
  temp = population_quantities(factors0)
  U_inv = temp[[1]]
  Psi = temp[[2]]
  Q = temp[[3]]
  H = temp[[4]]
  sigma2 = temp[[5]]
  sigma4 = temp[[6]]
  Sigma_lambda = temp[[7]]
  
  omega_hat = (t(factors0) %*% M_1T %*% factors0)/TT
  omega_hat_inv = solve(omega_hat)
  
  Ue = diag(sigma4, TT^2)
  for (j in 1:TT)
    for (k in 1:TT)
      Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
  
  D1 = (IvecT %x% Imatk) %*% omega_hat_inv %*% (t(factors0) %*% (ImatT[, t] - IvecT))
  D2 = ((ImatT[, t] - IvecT) %x% Imatk) %*% omega_hat_inv %*% (t(factors0) %*% IvecT)
  D31 = (1/TT * M_1T %*% factors0 %*% omega_hat_inv %*% t(factors0)) %x% (omega_hat_inv %*% t(factors0))
  D32 = (IvecT %x% (ImatT[, t] - IvecT)) + ((ImatT[, t] - IvecT) %x% IvecT)
  D = (ImatT %x% solve(H)) %*% (1/rf[t]*(- D1 - D2 + D31 %*% D32)) 
  
  E1 = (ImatT %x% (U_inv %*% (t(factors0 %*% H)/TT))) %*% Ue %*% (ImatT %x% ((factors0 %*% H/TT) %*% U_inv))
  E2 = sigma2*ImatT %x% U_inv
  E3 = (factors0 %*% Sigma_lambda %*% t(factors0)) %x% (sigma2/TT*U_inv^2)
  E4 = ((sigma2/TT* factors0 %*% H %*% U_inv) %x% t(factors0 %*% H)) %*% K_rT
  E5 = K_Tr %*% ((sigma2/TT*U_inv %*% t(factors0 %*% H)) %x% (factors0 %*% H))
  E = E1 + E2 + E3 + E4 + E5
  
  out = (t(D) %*% E %*% D)
  
  return(out)
}
SDFtilde_acm_fun = function(t, k, factors, eigenvalue, valuef, sigma4)
{
  TT = nrow(factors)

  rf = rep(1, TT)
  
  IvecT = matrix(1, ncol = 1, nrow = TT)/TT
  ImatT = diag(1, TT)
  Imatk = diag(1, k)
  M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
  K_rT = commutation_matrix(k, TT)
  K_Tr = commutation_matrix(TT, k)
  
  gamma = apply(factors, 2, mean)
  omega = diag(1, k) - gamma %*% t(gamma)
  omega_inv = solve(omega)
  sigma2 = TT/(TT - k)*valuef
  U = diag(eigenvalue, k) - sigma2/TT*Imatk
  U_inv = solve(U)
  sigma_lambda = U
  Ue = diag(sigma4, TT^2)
  for (j in 1:TT)
    for (k in 1:TT)
      Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
  
  D1 = (IvecT %x% Imatk) %*% omega_inv %*% (t(factors) %*% (ImatT[, t] - IvecT))
  D2 = ((ImatT[, t] - IvecT) %x% Imatk) %*% omega_inv %*% (t(factors) %*% IvecT)
  D31 = (1/TT * M_1T %*% factors %*% omega_inv %*% t(factors)) %x% (omega_inv %*% t(factors))
  D32 = (IvecT %x% (ImatT[, t] - IvecT)) + ((ImatT[, t] - IvecT) %x% IvecT)
  D = 1/rf[t]*(- D1 - D2 + D31 %*% D32)
  
  E1 = (ImatT %x% (U_inv %*% (t(factors)/TT))) %*% Ue %*% (ImatT %x% ((factors/TT) %*% U_inv))
  E2 = sigma2*ImatT %x% U_inv
  E3 = (factors %*% sigma_lambda %*% t(factors)) %x% (sigma2/TT*U_inv^2)
  E4 = ((sigma2/TT*factors %*% U_inv) %x% t(factors)) %*% K_rT
  E5 = K_Tr %*% ((sigma2/TT*U_inv %*% t(factors)) %x% factors)
  E = E1 + E2 + E3 + E4 + E5
  
  out = (t(D) %*% E %*% D)
  
  return(out)
}

clusterExport(cl, c("OLS", "nfac_selectioncriteria", "nfac_optimization", "extraction", "sqrtm",
                    "Ftilde_acm_fun", "RiskPremiatilde_fun", "RiskPremiatilde_acm_fun", "SDFtilde_fun", "SDFtilde_acm_fun",
                    "commutation_matrix", "population_quantities", "eigenvector_assumption","Ftilde_acm_funbai"))

Montecarlo = function(i, factors0, N, w = c("factors", "riskpremia", "sdf"))
{
  nfac0 = ncol(factors0)
  TT = nrow(factors0)
  
  # generazione dati
  loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  e = matrix(rnorm(TT*N, 0, 1), nrow = TT, ncol = N)
  
#  loadings = loadings - rep(1, N) %*% t(apply(loadings, 2, mean))     # delete these line at the end
#  loadings = loadings %*% solve(chol(cov(loadings) *(N - 1)/ N))      # delete these line at the end
#  t(loadings) %*% loadings / N                                        # delete these line at the end
  
  common = factors0 %*% t(loadings)
  X = common + sqrt(nfac0)*e
  
  #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
  nfac_hat = nfac0
  
  temp = nfac_optimization(data = X, kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  valuef = temp[[1]]
  eigenv = temp[[2]]
  loadmean = temp[[3]]
  loadvar = temp[[4]]
  factors = as.matrix(temp[[5]])
  sigma4 = temp[[6]]
  gamma2pass = temp[[7]]
  rm(temp)
  
  valuef = valuef[nfac_hat]
  eigenv = eigenv[1:nfac_hat]
  factors = as.matrix(factors[,1:nfac_hat])
  sigma4 = sigma4[nfac_hat]
  
  sigma2 = TT/(TT - nfac_hat)*valuef
  Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
  Utilde_inv = solve(Utilde)
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  temp = population_quantities(factors0)
  Psi = temp[[2]]

  Ltilde = rep(1, nfac0)
  for (i in 1:nfac0) 
  {
    t = 1
    while ((sign(round(Psi[t, i], 10)) == 0) | (sign(round(Psitilde[t,i], 10)) == 0)) t = t + 1
    if (sign(Psi[t,i]) != sign(Psitilde[t,i])) Ltilde[i] = -Ltilde[i]
  }
  
  factors = as.matrix(factors %*% diag(Ltilde, nfac0))
  Qtilde = t((1/TT * t(factors0) %*% factors))
  Htilde = (1/N * t(loadings) %*% loadings)  %*% (1/TT * t(factors0) %*% factors) %*% Utilde_inv
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  # factors
  Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT))
  if (any(w == "factors"))
  {
    for (t in 1:TT) Ftilde_acm[, , t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  }
  # risk premia
  RiskPremiatilde = rep(NA, nfac0)
  RiskPremiatilde_acm = matrix(NA, nfac0, nfac0)
  if (any(w == "riskpremia"))
  {
    RiskPremiatilde = RiskPremiatilde_fun(nfac_hat, factors)
    RiskPremiatilde_acm = RiskPremiatilde_acm_fun(nfac_hat, factors, eigenv, valuef, sigma4)
  }
  # SDF
  SDFtilde = rep(NA, TT)
  SDFtilde_acm = rep(NA, TT)
  if (any(w == "sdf"))
  {
    SDFtilde = SDFtilde_fun(nfac_hat, factors)
    for (t in 1:TT) SDFtilde_acm[t] = SDFtilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  }
  
  out = list(Qtilde, Htilde, Psitilde, Utilde_inv, factors, Ftilde_acm, RiskPremiatilde, RiskPremiatilde_acm, SDFtilde, SDFtilde_acm)
  return(out)
}

# time varying loadings: random walk

Montecarlo2 = function(i, seed0, factors0, hNT, N, w = c("factors", "riskpremia", "sdf"))
{
  #set.seed(seed0)
  nfac0 = ncol(factors0)
  TT = nrow(factors0)
  
  # generazione dati
  sigmaz=1
  loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  t=2
  xi=matrix(0,(N*nfac0),TT)
  for (t in 2:TT) 
  {
    
    xi[,t]=xi[,t-1]+rnorm((N*nfac0), 0, sigmaz)
   
      
  }
  
  
  e = matrix(rnorm(TT*N, 0, 1), nrow = TT, ncol = N)
  
  #  loadings = loadings - rep(1, N) %*% t(apply(loadings, 2, mean))     # delete these line at the end
  #  loadings = loadings %*% solve(chol(cov(loadings) *(N - 1)/ N))      # delete these line at the end
  #  t(loadings) %*% loadings / N                                        # delete these line at the end
  
  #common = factors0 %*% t(loadings)
  
  #hNT=1/log(log(N))
  hNT=hNT
  
  common=matrix(0,TT,N)
  t=1
  for (t in 1:TT)
  {
    common[t,]=(factors0[t,])%*%
      (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)))


  }
  
  X = common + sqrt(nfac0)*e
  
  #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
  nfac_hat = nfac0
  
  temp = nfac_optimization(data = X, kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  valuef = temp[[1]]
  eigenv = temp[[2]]
  loadmean = temp[[3]]
  loadvar = temp[[4]]
  factors = as.matrix(temp[[5]])
  
  sigma4 = temp[[6]]
  gamma2pass = temp[[7]]
  loadingshat=t(X)%*%factors/TT
  rm(temp)
  
  valuef = valuef[nfac_hat]
  eigenv = eigenv[1:nfac_hat]
  factors = as.matrix(factors[,1:nfac_hat])
  loadingshat = as.matrix(loadingshat[,1:nfac_hat])
  loadings = as.matrix(loadings[,1:nfac_hat])
  
  sigma4 = sigma4[nfac_hat]
  
  sigma2 = TT/(TT - nfac_hat)*valuef
  Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
  Utilde_inv = solve(Utilde)
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  temp = population_quantities(factors0)
  Psi = temp[[2]]
  
  Ltilde = rep(1, nfac0)
  for (i in 1:nfac0) 
  {
    t = 1
    while ((sign(round(Psi[t, i], 10)) == 0) | (sign(round(Psitilde[t,i], 10)) == 0)) t = t + 1
    if (sign(Psi[t,i]) != sign(Psitilde[t,i])) Ltilde[i] = -Ltilde[i]
  }
  
  factors = as.matrix(factors %*% diag(Ltilde, nfac0))
  loadings = as.matrix(loadings %*% diag(Ltilde, nfac0))
  
  Qtilde = t((1/TT * t(factors0) %*% factors))
  Htilde = (1/N * t(loadings) %*% loadings)  %*% (1/TT * t(factors0) %*% factors) %*% Utilde_inv
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  # factors
  Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT))
  Ftilde_acmbai = array(NA, dim = c(nfac0, nfac0, TT))
  
  if (any(w == "factors"))
  {
    for (t in 1:TT) Ftilde_acm[, , t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
    for (t in 1:TT) Ftilde_acmbai[, , t] = Ftilde_acm_funbai(t, nfac_hat, factors, eigenv, valuef, sigma4)
    
    }
  # risk premia
  RiskPremiatilde = rep(NA, nfac0)
  RiskPremiatilde_acm = matrix(NA, nfac0, nfac0)
  if (any(w == "riskpremia"))
  {
    RiskPremiatilde = RiskPremiatilde_fun(nfac_hat, factors)
    RiskPremiatilde_acm = RiskPremiatilde_acm_fun(nfac_hat, factors, eigenv, valuef, sigma4)
  }
  # SDF
  SDFtilde = rep(NA, TT)
  SDFtilde_acm = rep(NA, TT)
  if (any(w == "sdf"))
  {
    SDFtilde = SDFtilde_fun(nfac_hat, factors)
    for (t in 1:TT) SDFtilde_acm[t] = SDFtilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  }
  
  out = list(Qtilde, Htilde, Psitilde, Utilde_inv, factors, Ftilde_acm, RiskPremiatilde, RiskPremiatilde_acm, SDFtilde, SDFtilde_acm,loadingshat,loadings,Ftilde_acmbai,loadmean)
  return(out)
}

# IPCA with two instruments orthogonal plus local-PCA estimator

 Montecarlo3 = function(i, seed0 ,factors0, hNT, N, w = c("factors", "riskpremia", "sdf"))
{
  
  #set.seed(seed0)
  nfac0 = ncol(factors0)
  TT = nrow(factors0)
  
  
  sigmaz=1
  loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  t=2
  xi=matrix(0,(N*nfac0),TT)
  for (t in 2:TT)
  {

    xi[,t]=xi[,t-1]+rnorm((N*nfac0), 0, sigmaz)


  }
  
  
  
  e = matrix(rnorm(TT*N, 0, 1), nrow = TT, ncol = N)
  
  #  loadings = loadings - rep(1, N) %*% t(apply(loadings, 2, mean))     # delete these line at the end
  #  loadings = loadings %*% solve(chol(cov(loadings) *(N - 1)/ N))      # delete these line at the end
  #  t(loadings) %*% loadings / N                                        # delete these line at the end
  
  #common = factors0 %*% t(loadings)
  
  hNT=hNT #1/log(log(N))
  #hNT=10
  
  # generazione dati
  #  number char = 2
  nfacchar=2
  TTT=10*TT
  t=2
  meanchar=rnorm(nfacchar,0,1)#zero?
  char=matrix(0,(N*nfacchar),TTT)
  for (t in 2:TTT) 
  {
    char[,t]=(1-0.8)*(meanchar%x%matrix(1,N,1))  +  0.8*char[,t-1]+rnorm(N*nfacchar,0,1)
  }  
  
  gammab1=matrix(rnorm((nfac0*nfacchar),0,1),nfacchar,nfac0)
  #gammab0=rnorm((nfac0*nfacchar),0,1)
  
  charfinal=char[,(TTT-TT+1):TTT]
  
  # ortho char
  orthchar=array(NA, dim = c(nfacchar, N , TTT))
    t=2
  for (t in 2:TTT)
  {
    mchar=matrix(char[,t],nfacchar,N)
    orthchar[,,t]=solve(sqrtm(mchar%*%t(mchar)))%*%mchar
  }
  
    #   orthchar[,,t]%*%t(orthchar[,,t])
  
  common=matrix(0,TT,N)
  loadingsIPCA=array(0,dim=c(N,nfac0,TT))
  loadingsIPCA_wrong=array(0,dim=c(N,nfac0,TT))
  
  # common[1,]=t(factors0[1,])%*%
  #   (t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #      hNT*t(gammab1)%*%matrix(char[,(TTT-TT)],nfacchar,N))    
  
  common[1,]=t(factors0[1,])%*%
    (t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
       hNT*t(gammab1)%*%orthchar[,,(TTT-TT)])
  
  loadingsIPCA[,,1]=as.matrix(t(t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
      hNT*t(gammab1)%*%orthchar[,,(TTT-TT)]),N,nfac0)
  
  
  t=2
  for (t in 2:TT)
  {
    
    # common[t,]=(factors0[t,])%*%
    #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
    #      hNT*t(gammab1)%*%matrix(charfinal[,t-1],nfacchar,N))    
    # 
    
    common[t,]=t(factors0[t,])%*%
      (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
         hNT*t(gammab1)%*%orthchar[,,(TTT-TT+t-1)])
    
    loadingsIPCA[,,t]=t(    (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
        hNT*t(gammab1)%*%orthchar[,,(TTT-TT+t-1)]) )
    
    
  }
  
  X = common + sqrt(nfac0)*e
  X=as.matrix(X,TT,N)
  
  #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
 
  #  estimates IPCA
    USU=matrix(0,nfacchar,nfacchar)
    tt=1
    for (tt in 1:TT)
    {
      USU=USU +   t(X[tt,]%*%t(orthchar[,,(TTT-TT+tt-1)]))%*%X[tt,]%*%t(orthchar[,,(TTT-TT+tt-1)])
    }
    
  USU.svd=svd(USU)
  
  
  aaa= as.matrix(USU.svd$u,nfacchar,nfacchar)
  gammab1tilde=aaa[,1:nfac0] #  kc times k
  gammab1tilde_wrong=aaa[1,1:nfac0] #  kcwrong=1 times k
  
  
  FtildeIPCA=matrix(0,TT,nfac0)
  FtildeIPCA_wrong=matrix(0,TT,nfac0)
  
  tt=1
  for (tt in 1:TT)
  {
    FtildeIPCA[tt,]=t(X[tt,]%*%t(orthchar[,,(TTT-TT+tt-1)])%*%gammab1tilde) 
    FtildeIPCA_wrong[tt,]=t(X[tt,]%*%t(matrix(orthchar[1,,(TTT-TT+tt-1)],1,N))%*%gammab1tilde_wrong) 
    
    
    # commontildeIPCA[tt,]=t(FtildeIPCA[tt,])%*%
    #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
    #      hNT*t(gammab1)%*%orthchar[,,(TTT-TT+t-1)])
    # 
    
    #t(gammab1tilde)%*%orthchar[,,(TTT-TT+tt-1)]%*%t(X[tt,])
  }

  loadingstildeIPCA=array(0,dim=c(N,nfac0,T))
  loadingstildeIPCA_wrong=array(0,dim=c(N,nfac0,T))
  
  
  t=1
  for (i in t:TT)
  {
   loadingstildeIPCA[,,t]=t(orthchar[,,(TTT-TT+t-1)])%*%gammab1tilde 
   loadingstildeIPCA_wrong[,,t]=as.matrix(orthchar[1,,(TTT-TT+t-1)],N,1)*gammab1tilde_wrong
     
    #t(gammab1tilde)%*%orthchar[,,(TTT-TT+tt-1)]%*%t(X[tt,])
  }
  
  meanloadings0=apply(loadingsIPCA,c(1,2),mean)
# this is a N x nfac0 - averaged across t=1,...,TT=5  
  
  nfac_hat = nfac0
  
  temp = nfac_optimization(data = X, kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  valuef = temp[[1]]
  eigenv = temp[[2]]
  loadmean = temp[[3]]
  loadvar = temp[[4]]
  factors = as.matrix(temp[[5]])
  
  sigma4 = temp[[6]]
  gamma2pass = temp[[7]]
  loadingshat=t(X)%*%factors/TT
  rm(temp)
  
  valuef = valuef[nfac_hat]
  eigenv = eigenv[1:nfac_hat]
  factors = as.matrix(factors[,1:nfac_hat])
  loadingshat = as.matrix(loadingshat[,1:nfac_hat])
  #loadings = as.matrix(loadings[,1:nfac_hat])
  
  sigma4 = sigma4[nfac_hat]
  
  sigma2 = TT/(TT - nfac_hat)*valuef
  Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
  Utilde_inv = solve(Utilde)
#  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(meanloadings0) %*% meanloadings0)))
  
  temp = population_quantities(factors0)
  Psi = temp[[2]]
  
  Ltilde = rep(1, nfac0)
  for (i in 1:nfac0) 
  {
    t = 1
    while ((sign(round(Psi[t, i], 10)) == 0) | (sign(round(Psitilde[t,i], 10)) == 0)) t = t + 1
    if (sign(Psi[t,i]) != sign(Psitilde[t,i])) Ltilde[i] = -Ltilde[i]
  }
  
  factors = as.matrix(factors %*% diag(Ltilde, nfac0))
  loadings=meanloadings0 # important! approx equal to constant portion of loadings
  loadings = as.matrix(loadings %*% diag(Ltilde, nfac0))
  
  
  Qtilde = t((1/TT * t(factors0) %*% factors))
  Htilde = (1/N * t(loadings) %*% loadings)  %*% (1/TT * t(factors0) %*% factors) %*% Utilde_inv
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  # factors
  Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT))
  Ftilde_acmbai = array(NA, dim = c(nfac0, nfac0, TT))
  
  if (any(w == "factors"))
  {
    for (t in 1:TT) Ftilde_acm[, , t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
    for (t in 1:TT) Ftilde_acmbai[, , t] = Ftilde_acm_funbai(t, nfac_hat, factors, eigenv, valuef, sigma4)
    
  }
  # risk premia
  RiskPremiatilde = rep(NA, nfac0)
  RiskPremiatilde_acm = matrix(NA, nfac0, nfac0)
  if (any(w == "riskpremia"))
  {
    RiskPremiatilde = RiskPremiatilde_fun(nfac_hat, factors)
    RiskPremiatilde_acm = RiskPremiatilde_acm_fun(nfac_hat, factors, eigenv, valuef, sigma4)
  }
  # SDF
  SDFtilde = rep(NA, TT)
  SDFtilde_acm = rep(NA, TT)
  if (any(w == "sdf"))
  {
    SDFtilde = SDFtilde_fun(nfac_hat, factors)
    for (t in 1:TT) SDFtilde_acm[t] = SDFtilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  }
  
  out = list(Qtilde, Htilde, Psitilde, Utilde_inv, factors, Ftilde_acm, RiskPremiatilde, RiskPremiatilde_acm, SDFtilde, SDFtilde_acm,loadingshat,
             loadings,Ftilde_acmbai,loadmean,FtildeIPCA,FtildeIPCA_wrong,
             loadingstildeIPCA,loadingstildeIPCA_wrong,loadingsIPCA,loadingsIPCA_wrong)
  
             
  return(out)
}
}





Montecarlo4 = function(i, q , eps1, indIPCA1, indIPCA2, nfacchar, seed0 ,factors0, TT, hNT, N, w = c("factors", "riskpremia", "sdf"))
{
  
  # this part including estimation IPCA does not roll
  
  #set.seed(seed0)
  eps1=eps1
  nfac0 = ncol(factors0) # long time series TTT
  TTT = nrow(factors0)
  TT=TT
  nfacchar=nfacchar
  indIPCA1=indIPCA1
  indIPCA2=indIPCA2
  
  
  
  sigmaz=1
  loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  t=2
  xi=matrix(0,(N*nfac0),TTT)
  for (t in 2:TTT)
  {
    
    xi[,t]=xi[,t-1]+rnorm((N*nfac0), 0, sigmaz)
    
    
  }
  
  #  assume  e_it = delta_i u_t + eta_it 
  # eta_it iid across time and i
  #  u_t iid across time
  # delta_i = N^eps1 for i=1....q and  delta_i = 1/N^(0.5+eps) i=q+1,...,N
  # q binomial or exponential or uniform or constant?
  
  q=q#N/10
  #eps1=0.25/2 # should break when > 1/4
  delta=matrix(1,N,1)
  delta[1:q,]=delta[1:q]*N^eps1
  delta[(q+1):N,]=delta[(q+1):N,]*N^(-(0.5+eps1))
  e=rnorm(TTT,0,1)%*%t(delta) + matrix(rnorm(TTT*N, 0, 1), nrow = TTT, ncol = N)
  e=e/sqrt(mean(diag(var(e))))
  
  
  
  
#  e = matrix(rnorm(TTT*N, 0, 1), nrow = TTT, ncol = N)
  
  #  loadings = loadings - rep(1, N) %*% t(apply(loadings, 2, mean))     # delete these line at the end
  #  loadings = loadings %*% solve(chol(cov(loadings) *(N - 1)/ N))      # delete these line at the end
  #  t(loadings) %*% loadings / N                                        # delete these line at the end
  
  #common = factors0 %*% t(loadings)
  
  hNT=hNT #1/log(log(N))
  #hNT=10
  
  # generazione dati
  #  number char = 2
  #nfacchar=10
  TTT1=2*TTT
  t=2
  meanchar=rnorm(nfacchar,0,1)#zero?
  char=matrix(0,(N*nfacchar),TTT1)
  for (t in 2:TTT1) 
  {
    char[,t]=(1-0.8)*(meanchar%x%matrix(1,N,1))  +  0.8*char[,t-1]+rnorm(N*nfacchar,0,1)
  }  
  
  gammab1=matrix(rnorm((nfac0*nfacchar),0,1),nfacchar,nfac0)
  gammab1=gammab1%*%solve(sqrtm(t(gammab1)%*%gammab1))
  #gammab0=rnorm((nfac0*nfacchar),0,1)
  
  charfinal=char[,(TTT1-TTT+1):TTT1]
  
  # ortho char
  orthchar=array(NA, dim = c(nfacchar, N , TTT1))
  t=2
  for (t in 2:TTT1)
  {
    mchar=matrix(char[,t],nfacchar,N)
    orthchar[,,t]=solve(sqrtm(mchar%*%t(mchar)))%*%mchar
  }
  
  #   orthchar[,,t]%*%t(orthchar[,,t])
  
  common=matrix(0,TTT,N)
  loadingsIPCA=array(0,dim=c(N,nfac0,TTT))
  loadingsIPCA_wrong=array(0,dim=c(N,nfac0,TTT))
  loadingsIPCA2=array(0,dim=c(N,nfac0,TTT))
  
  # common[1,]=t(factors0[1,])%*%
  #   (t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #      hNT*t(gammab1)%*%matrix(char[,(TTT-TT)],nfacchar,N))    
  
    # 
    # (0*t(loadings + 0*hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
    #    hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)])
    # 
  loadingsIPCA[,,1]=
    as.matrix(t(t(indIPCA1*loadings + indIPCA2*hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
                                  hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)]),N,nfac0)
  
  loadingsIPCA2[,,1]=as.matrix(t(hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)]),N,nfac0)
  
  common[1,]=t(factors0[1,])%*%t(loadingsIPCA[,,1])
    
  t=2
  for (t in 2:TTT)
  {
    
    # common[t,]=(factors0[t,])%*%
    #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
    #      hNT*t(gammab1)%*%matrix(charfinal[,t-1],nfacchar,N))    
    # 
    
      
      # (t(0*loadings + 0*hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
      #    hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)])
    
    loadingsIPCA[,,t]=t(    (t(indIPCA1*loadings + indIPCA2*hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
                               hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)]) )
    
    loadingsIPCA2[,,t]=t(      hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)] )
  
    common[t,]=t(factors0[t,])%*%t(loadingsIPCA[,,t])
    
      
  }
  
  ss=mean(diag(var(common)))
  
  #X = common + sqrt(nfac0)*e*sqrt(ss)
  X = common + e*sqrt(ss)*2.5
  X=matrix(X,TTT,N)
  # QQQ
  #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
  
  #  estimates IPCA
  USU=matrix(0,nfacchar,nfacchar)
  tt=1
  for (tt in 1:TTT)
  {
    USU=USU +   t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)]))%*%X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])
  }
  
  USU.svd=svd(USU)
  
  
  
  aaa= as.matrix(USU.svd$u,nfacchar,nfacchar)
  gammab1tilde=matrix(aaa[,1:nfac0],nfacchar,nfac0) #  kc times k
  gammab1tilde_wrong=matrix(aaa[1,1:nfac0],1,nfac0) #  kcwrong=1 times k
  
  
  #Ltilde = rep(1, nfac0)
  # for (i in 1:nfac0) 
  # {
  #   t = 1
  #   while ((sign(round(gammab1[t, i], 10)) == 0) | (sign(round(gammab1tilde[t,i], 10)) == 0)) t = t + 1
  #   if (sign(gammab1[t,i]) != sign(gammab1tilde[t,i])) Ltilde[i] = -Ltilde[i]
  # }
  # 

  Ltilde=matrix(0,nfac0,nfac0)

for (i in 1:nfac0)
{
  if (gammab1tilde[1,i]*gammab1[1,i]>0) { Ltilde[i,i]=1} else {Ltilde[i,i]=-1}
}



  gammab1tilde=gammab1tilde%*%Ltilde
  gammab1tilde_wrong=gammab1tilde_wrong%*%Ltilde
  
   
  #HtildeIPCA=svd(t(gammab1)%*%gammab1)$u
  #gammmab1tilderot=gammab1tilde%*%t(HtildeIPCA)
  
  FtildeIPCA=matrix(0,TTT,nfac0)
  FtildeIPCA_wrong=matrix(0,TTT,nfac0)
  
  tt=1
  for (tt in 1:TTT)
  {
    FtildeIPCA[tt,]=t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])%*%gammab1tilde) 
    FtildeIPCA_wrong[tt,]=t(X[tt,]%*%t(matrix(orthchar[1,,(TTT1-TTT+tt-1)],1,N))%*%gammab1tilde_wrong) 
    
    
    # commontildeIPCA[tt,]=t(FtildeIPCA[tt,])%*%
    #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
    #      hNT*t(gammab1)%*%orthchar[,,(TTT-TT+t-1)])
    # 
    
    #t(gammab1tilde)%*%orthchar[,,(TTT-TT+tt-1)]%*%t(X[tt,])
  }
  
  loadingstildeIPCA=array(0,dim=c(N,nfac0,TTT))
  loadingstildeIPCA_wrong=array(0,dim=c(N,nfac0,TTT))
  commontildeIPCA=array(0,dim=c(TTT,N))
  commontildeIPCAwrong=array(0,dim=c(TTT,N))
  
    # remember that in IPCA one only estimates a bit of the
  # total loadings, that is the one that is related to gammab1
  
  t=1
  for (t in 1:TTT)
  {
    loadingstildeIPCA[,,t]=t(orthchar[,,(TTT1-TTT+t-1)])%*%gammab1tilde 
    loadingstildeIPCA_wrong[,,t]=as.matrix(orthchar[1,,(TTT1-TTT+t-1)],N,1)%*%gammab1tilde_wrong
    
    commontildeIPCA[t,]=FtildeIPCA[t,]%*%t(loadingstildeIPCA[,,t])
    commontildeIPCAwrong[t,]=FtildeIPCA[t,]%*%t(loadingstildeIPCA_wrong[,,t])
    
  }
  
  
  
  UtildeIPCA=diag(USU.svd$d[1:nfac0],nfac0)
  UtildeIPCA_inv=solve(UtildeIPCA)
  
  HtildeIPCA=matrix(0,nfac0,nfac0)
  HtildeIPCA = ( t(gammab1) %*% gammab1)  %*%(1/TTT * t(factors0) %*%FtildeIPCA) %*%  UtildeIPCA_inv
  
  
  
  # # OOS  IPCA
  USU=matrix(0,nfacchar,nfacchar)
  gammab1tildeo=array(0,dim=c(nfacchar,nfac0,TTT))
  gammab1tildeowrong=array(0,dim=c(1,nfac0,TTT))
  

  tt=1
  for (tt in 1:(TT-1))
  {
    USU=USU +   t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)]))%*%X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])
  }

  aaa= as.matrix(USU.svd$u,nfacchar,nfacchar)
  gammab1tildeo[,,TT-1]=matrix(aaa[,1:nfac0],nfacchar,nfac0) #  kc times k
  gammab1tildeowrong[,,TT-1]=matrix(aaa[1,1:nfac0],nfac0)#  kc times k
  

  FtildeIPCAo=matrix(NA,TTT,nfac0)
  FtildeIPCAowrong=matrix(NA,TTT,nfac0)
  loadingstildeIPCAo=array(NA,dim=c(N,nfac0,TTT))
  loadingstildeIPCAowrong=array(NA,dim=c(N,nfac0,TTT))
  commontildeIPCAo=matrix(NA,TTT,N)
  commontildeIPCAowrong=matrix(NA,TTT,N)
  


  tt=TT
  for (tt in TT:TTT)
  {
    USU=USU +   t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)]))%*%X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])

    USU.svd=svd(USU)

    aaa= as.matrix(USU.svd$u,nfacchar,nfacchar)
    gammab1tildeo[,,tt]=matrix(aaa[,1:nfac0],nfacchar,nfac0) #  kc times k
    gammab1tildeowrong[,,tt]=matrix(aaa[1,1:nfac0],1,nfac0) #  kc times k

    FtildeIPCAo[tt,]=t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])%*%gammab1tildeo[,,tt-1])
    FtildeIPCAowrong[tt,]=t(X[tt,]%*%t(matrix(orthchar[1,,(TTT1-TTT+tt-1)],1,N))%*%gammab1tildeowrong[,,tt-1])
    
    loadingstildeIPCAo[,,tt]=t(orthchar[,,(TTT1-TTT+tt-1)])%*%gammab1tildeo[,,tt-1]
    loadingstildeIPCAowrong[,,tt]=matrix(orthchar[1,,(TTT1-TTT+tt-1)],N,1)%*%matrix(gammab1tildeowrong[,,tt-1],1,nfac0)

   commontildeIPCAo[tt,]=matrix(FtildeIPCAo[tt,],1,nfac0)%*%t(matrix(loadingstildeIPCAo[,,tt],N,nfac0))
    commontildeIPCAowrong[tt,]=matrix(FtildeIPCAowrong[tt,],1,nfac0)%*%t(loadingstildeIPCAowrong[,,tt])


  }
  # 
  # # find HtildeIPCA? not need as gammab is fixed
  
  # rolling estimation of local PCA
  
  nfac_hat = nfac0  
  
  RO=TTT-TT+1
  #TT = nrow(X)
  #N = ncol(X)
  # eigenv = matrix(NA,nfac_hat,RO) #eigen(1/(TT*N) * X %*% t(X))$values[1:nk]
  # valuef = array(NA,dim=c(nfac_hat,RO)) #rep(NA, nk)
  #                m_loadings =matrix(0, ncol = nk, nrow = nk)
  #                s_loadings =matrix(NA, ncol = nk, nrow = nk)
                  factors = array(NA,dim=c(TT,nfac0,TTT))#matrix(NA, ncol = nk, nrow = TT)
                 
                  # sigma_4 = matrix(NA,nfac_hat,RO)#rep(NA, nk)
                 # gamma_twopass = array(NA,dim=c(nfac_hat,nfac_hat,RO)) #matrix(NA, ncol = nk, nrow = nk)
                  loadingshat = array(NA,dim=c(N,nfac0,TTT))#matrix(0, ncol = nk, nrow = nk)
                  factorshat = array(NA,dim=c(TTT,nfac0))#matrix(0, ncol = nk, nrow = nk)
                  Htilde=array(NA,dim=c(nfac0,nfac0,TTT))
               
                  loadingshato = array(NA,dim=c(N,nfac0,TTT))#matrix(0, ncol = nk, nrow = nk)
                  factorshato = array(NA,dim=c(TTT,nfac0))#matrix(0, ncol = nk, nrow = nk)
                  
                  
  
  roll=TT
  for (roll in TT:TTT)
  {
    
    
    
    #temp = nfac_optimization(data = X[(roll-TT+1):roll,], kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
    
  #  TT = nrow(X)
  #  N = ncol(X)
    nk=nfac_hat
    eig = eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$values[1:nk]
    
    val = rep(NA, nk)
    m_loadings = matrix(0, ncol = nk, nrow = nk)
    s_loadings = matrix(NA, ncol = nk, nrow = nk)
    #factors = matrix(NA, ncol = nk, nrow = TT)
    sigma_4 = rep(NA, nk)
    gamma_twopass = matrix(NA, ncol = nk, nrow = nk)

    k=1
    for (k in 1:nk)
    {
      temp = sqrt(TT) * as.matrix(eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$vectors[,1:k])
      loadings = 1/TT *(t(X[(roll-TT+1):roll,]) %*% temp)
      m_loadings[1:k, k] = apply(loadings, 2, mean)
      val[k] = 1/(N*TT) * sum((X[(roll-TT+1):roll,] - temp %*% t(loadings))^2)
      temp1 = sum(diag((temp %*% t(temp))/TT)^2)
      temp1 = 3 + 27/TT*temp1 + 18*k/TT
      temp2 = 1/(N*TT) * sum((X[(roll-TT+1):roll,] - temp %*% t(loadings))^4)
      sigma_4[k] = temp2/temp1
      gamma_twopass[1:k,k] = solve(t(loadings) %*% loadings) %*% t(loadings) %*% apply(X[(roll-TT+1):roll,], 2, mean)
    
    #factors[1:TT, 1:k] = temp
    s_loadings[1:k, 1:k] = (t(loadings) %*% loadings)/N
    }
    
    valuef = val[nfac_hat]
    eigenv = eig[1:nfac_hat]
    s_loadings=s_loadings[1:nfac_hat,1:nfac_hat]
    gamma_twopass=gamma_twopass[1:nfac_hat,1:nfac_hat]
    m_loadings=m_loadings[1:nfac_hat,1:nfac_hat]
    
    #temp = nfac_optimization(data = X[roll-TT+1:roll,], kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
    # valuef = temp[[1]]
    # eigenv = temp[[2]]
    # loadmean= temp[[3]]
    # loadvar = temp[[4]]
    #factors = as.matrix(temp[[5]])
    factors[,,roll] =  sqrt(TT) * as.matrix(eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$vectors[,1:nfac_hat])
    
    
    #sigma4 = temp[[6]]
    loadingshat[,,roll]=t(X[(roll-TT+1):roll,])%*%factors[,,roll]/TT
 
    
    roll1=roll-1
    loadingshato[,,roll]=  sqrt(N) * as.matrix(eigen(1/(TT*N) *  t(X[(roll1-TT+1):roll1,])%*%X[(roll1-TT+1):roll1,])$vectors[,1:nfac_hat])
    factorshato[roll,]=   t(matrix(loadingshato[,,roll],N,nfac0))%*%t(matrix(X[roll,],1,N))/N
    


    
       #gamma2pass= temp[[7]]
   # rm(temp)
    
     #valuef = valuef[nfac_hat]
     #eigenv = eigenv[1:nfac_hat]
   # factors = as.matrix(factors[,1:nfac_hat])
  #  loadingshat = as.matrix(loadingshat[,1:nfac_hat])
  #  loadings = as.matrix(loadings[,1:nfac_hat])
    
     sigma_4 = sigma_4[nfac_hat]
    # 
     sigma2 = TT/(TT - nfac_hat)*valuef
     Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
     Utilde_inv = solve(Utilde)
     Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(matrix(factors0[(roll-TT+1):roll,],TT,nfac0)) %*% factors[,,roll])) %*% sqrtm((1/N * t(loadings) %*% loadings)))
    # 
    temp = population_quantities(matrix(factors0[(roll-TT+1):roll,1:nk],TT,nk))
    Psi = temp[[2]]
    
    Ltilde = rep(1, nfac0)
    for (i in 1:nfac0) 
    {
      t = 1
      while ((sign(round(Psi[t, i], 10)) == 0) | (sign(round(Psitilde[t,i], 10)) == 0)) t = t + 1
      if (sign(Psi[t,i]) != sign(Psitilde[t,i])) Ltilde[i] = -Ltilde[i]
    }
    
    factors[,,roll] = as.matrix(factors[,,roll] %*% diag(Ltilde, nfac0))
    loadingshat[,,roll] = as.matrix(loadingshat[,,roll] %*% diag(Ltilde, nfac0))
    
    factorshat[roll,]=factors[TT,,roll] # this is the last estimate from rolling, then glued together
    
     meanloadings0=apply(array(loadingsIPCA[,,(roll-TT+1):roll],dim=c(N,nfac0,TT)),c(1,2),mean)
     loadings=meanloadings0
    # Qtilde[,,roll] = t((1/TT * t(factors0[(roll-TT+1):roll,]) %*% factors[,,roll]))
     Htilde[,,roll] = (1/N * t(loadings) %*% loadings)  %*% 
       matrix((1/TT * t(matrix(factors0[(roll-TT+1):roll,1:nfac0],TT,nfac0)) %*% matrix(factors[,,roll],TT,nfac0)),nfac0,nfac0) %*% Utilde_inv
    
     
     
     # Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
    # 
    # # factors
    # Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT))
    # Ftilde_acmbai = array(NA, dim = c(nfac0, nfac0, TT))
    # 
    # if (any(w == "factors"))
    # {
    #   for (t in 1:TT) Ftilde_acm[, , t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
    #   for (t in 1:TT) Ftilde_acmbai[, , t] = Ftilde_acm_funbai(t, nfac_hat, factors, eigenv, valuef, sigma4)
    #   
    # }
    # # risk premia
    # RiskPremiatilde = rep(NA, nfac0)
    # RiskPremiatilde_acm = matrix(NA, nfac0, nfac0)
    # if (any(w == "riskpremia"))
    # {
    #   RiskPremiatilde = RiskPremiatilde_fun(nfac_hat, factors)
    #   RiskPremiatilde_acm = RiskPremiatilde_acm_fun(nfac_hat, factors, eigenv, valuef, sigma4)
    # }
    # # SDF
    # SDFtilde = rep(NA, TT)
    # SDFtilde_acm = rep(NA, TT)
    # if (any(w == "sdf"))
    # {
    #   SDFtilde = SDFtilde_fun(nfac_hat, factors)
    #   for (t in 1:TT) SDFtilde_acm[t] = SDFtilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
    # }
  }
  
  factorshat[1:TT,]=factors[,,TT]
  
commontilde=matrix(0,TTT,N)
commontildeo=matrix(0,TTT,N)

t=1
for (t in 1:TTT)
{ 
  
  commontilde[t,]=factorshat[t,]%*%t(loadingshat[,,t])
  commontildeo[t,]=factorshato[t,]%*%t(loadingshato[,,t])
  

}

# pricing errors

pe=matrix(0,TTT,N)
pei=matrix(0,TTT,N)
peiw=matrix(0,TTT,N)

pe2=matrix(0,TTT,1)
pei2=matrix(0,TTT,1)
peiw2=matrix(0,TTT,1)


peo=matrix(0,TTT,N)
peio=matrix(0,TTT,N)
peiwo=matrix(0,TTT,N)

peo2=matrix(0,TTT,1)
peio2=matrix(0,TTT,1)
peiwo2=matrix(0,TTT,1)





TT2=2*TT
t=TT2
for (t in TT2:TTT)
{  commontilde[t,]=factorshat[t,]%*%t(loadingshat[,,t])




pe[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontilde[(t-TT+1):t,],2,mean)
pe2[t]=pe[t,]%*%t(pe[t,])/N

pei[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontildeIPCA[(t-TT+1):t,],2,mean)
pei2[t]=pei[t,]%*%t(pei[t,])/N

peiw[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontildeIPCAwrong[(t-TT+1):t,],2,mean)
peiw2[t]=peiw[t,]%*%t(peiw[t,])/N


peo[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontilde[(t-TT):(t-1),],2,mean)
peo2[t]=peo[t,]%*%t(peo[t,])/N


peio[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontildeIPCAo[(t-TT+1):t,],2,mean)
peio2[t]=peio[t,]%*%t(peio[t,])/N

peiwo[t,]=apply(X[(t-TT+1):t,],2,mean) - apply(commontildeIPCAowrong[(t-TT+1):t,],2,mean)
peiwo2[t]=peiwo[t,]%*%t(peiwo[t,])/N



  
  }


# R2 
r2=0
r2i=0
r2iw=0
r2p=0
r2ip=0
r2iwp=0
r2o=0
r2io=0
r2iwo=0
r2po=0
r2ipo=0
r2iwpo=0



den=0

meanFtilde=apply(factorshat,2,mean)
meanFtildeIPCA=apply(FtildeIPCA,2,mean)
meanFtildeIPCAwrong=apply(FtildeIPCA_wrong,2,mean)


t=2*TT
for (t in (2*TT):TTT)
{  

den=den+matrix(X[t,],1,N)%*%matrix(t(X[t,]),N,1)
r2=r2+matrix((X[t,]-commontilde[t,]),1,N)%*%matrix(t(X[t,]-commontilde[t,]),N,1)
r2i=r2i+matrix((X[t,]-commontildeIPCA[t,]),1,N)%*%matrix(t(X[t,]-commontildeIPCA[t,]),N,1)
r2iw=r2iw+matrix((X[t,]-commontildeIPCAwrong[t,]),1,N)%*%matrix(t(X[t,]-commontildeIPCAwrong[t,]),N,1)

# this is kelly measure 1
r2p=r2p+(X[t,]-meanFtilde%*%t(loadingshat[,,t]))%*%t(X[t,]-meanFtilde%*%t(loadingshat[,,t]))
r2ip=r2ip+(X[t,]-meanFtildeIPCA%*%t(loadingstildeIPCA[,,t]))%*%t(X[t,]-meanFtildeIPCA%*%t(loadingstildeIPCA[,,t]))
r2iwp=r2iwp+(X[t,]-meanFtildeIPCAwrong%*%t(loadingstildeIPCA_wrong[,,t]))%*%t(X[t,]-meanFtildeIPCAwrong%*%t(loadingstildeIPCA_wrong[,,t]))

# this is kelly measure 2
r2o=r2o+matrix((X[t,]-commontildeo[t,]),1,N)%*%matrix(t(X[t,]-commontildeo[t,]),N,1)
r2io=r2io+matrix((X[t,]-commontildeIPCAo[t,]),1,N)%*%matrix(t(X[t,]-commontildeIPCAo[t,]),N,1)
r2iwo=r2iwo+matrix((X[t,]-commontildeIPCAowrong[t,]),1,N)%*%matrix(t(X[t,]-commontildeIPCAowrong[t,]),N,1)


meanFtildeo=matrix(apply(matrix(factorshato[(t-TT+1):t,],TT,nfac0),2,mean),nfac0,1)
meanFtildeIPCAo=matrix(apply(matrix(FtildeIPCAo[TT:t,],(t-TT+1),nfac0),2,mean),nfac0,1)
meanFtildeIPCAowrong=matrix(apply(matrix(FtildeIPCAowrong[TT:t,],(t-TT+1),nfac0),2,mean),nfac0,1)


#  this is oos versio of kelly 1
r2po=r2po+(X[t,]-t(meanFtildeo)%*%t(loadingshato[,,t]))%*%t(X[t,]-t(meanFtildeo)%*%t(loadingshato[,,t]))
r2ipo=r2ipo+(X[t,]-t(meanFtildeIPCAo)%*%t(loadingstildeIPCAo[,,t]))%*%t(X[t,]-t(meanFtildeIPCAo)%*%t(loadingstildeIPCAo[,,t]))
r2iwpo=r2iwpo+(X[t,]-t(meanFtildeIPCAowrong)%*%t(loadingstildeIPCAowrong[,,t]))%*%t(X[t,]-t(meanFtildeIPCAowrong)%*%t(loadingstildeIPCA_wrong[,,t]))





}
  
 r2=1-r2/den
 #r2
 r2i=1-r2i/den
 #r2i
 r2iw=1-r2iw/den
 #r2iw
 r2p=1-r2p/den
 #r2p
 r2ip=1-r2ip/den
 #r2ip
 r2iwp=1-r2iwp/den
 #r2iwp
 
 r2o=1-r2o/den
 #r2p
 r2io=1-r2io/den
 #r2ip
 r2iwo=1-r2iwo/den
 #r2iwp
 
 
 r2po=1-r2po/den
 #r2p
 r2ipo=1-r2ipo/den
 #r2ip
 r2iwpo=1-r2iwpo/den
 #r2iwp
 

 
 
    # out = list(Qtilde, Htilde, Psitilde, Utilde_inv, factors, Ftilde_acm, RiskPremiatilde, RiskPremiatilde_acm, SDFtilde, SDFtilde_acm,loadingshat,
    #            loadings,Ftilde_acmbai,loadmean,FtildeIPCA,FtildeIPCA_wrong,
    #            loadingstildeIPCA,loadingstildeIPCA_wrong,loadingsIPCA,loadingsIPCA_wrong)
  out = list(loadingsIPCA,FtildeIPCA,FtildeIPCA_wrong,loadingstildeIPCA,loadingstildeIPCA_wrong,loadingshat,factorshat,gammab1,gammab1tilde,
gammab1tilde_wrong,Htilde,common,commontildeIPCA,commontilde,loadingsIPCA2,commontildeIPCAwrong,pe2,pei2,peiw2,
r2,r2i,r2iw,r2p,r2ip,r2iwp,r2o,r2io,r2iwo,X,HtildeIPCA,
FtildeIPCAo,FtildeIPCAowrong,loadingstildeIPCAo,
loadingstildeIPCAowrong,commontildeIPCAo,commontildeIPCAowrong,r2po,r2ipo,r2iwpo,peo2,peio2,peiwo2)
    return(out)
  }

Nind=rbind(15,100,1000)
Tind=rbind(5,25)
pca2=rbind(0,1)
# Input
jj=1
for ( jj in 1:3)
{
  print(jj)
ii=1
for ( ii in 1:2)
{
  print(ii)
bb=1
for (bb in 1:2)
  {
  print(bb)
#  jj=1
#  ii=2
#  bb=1
cbind(jj,ii,bb)
N = Nind[jj]#10
TTT=100
TT=Tind[ii]#5
indIPCA1=0#pca1[aa]  # when 1 loadings with constant; when 0 only IPCA part
indIPCA2=pca2[bb]  # when 1 loadings with rw; when 0 only IPCA part
q=N/10#N-1 # between 1 and N-1.
eps1=0.2/2 # should break when > 1/4
seed0=123
hNT=1/(sqrt(N)*log(log(N)))
nfac0 = 3
nfacchar=10
meanF0 = rep(0, nfac0)  
SigmaF0 = diag((1:nfac0), nfac0)
nsimul = 100

# Application starts
set.seed(2)
F0 = matrix(NA, nrow = TTT, ncol = nfac0)
for (t in 1:TTT) F0[t,] = mvrnorm(1, meanF0, SigmaF0)
factors0=F0%*%solve(sqrtm(t(F0)%*%F0))

# inizializzazioni
{
Ftilde = array(NA, dim = c(TTT, nfac0, nsimul))
loadingstilde = array(NA, dim = c(N, nfac0, TTT, nsimul))
FtildeIPCA = array(NA, dim = c(TTT,nfac0,nsimul))
FtildeIPCAwrong = array(NA, dim = c(TTT,nfac0,nsimul))
loadingstildeIPCA= array(NA, dim = c(N,nfac0,TTT,nsimul))
loadingstildeIPCA_wrong=array(NA, dim = c(N,nfac0,TTT,nsimul))
loadingsIPCA= array(NA, dim = c(N,nfac0,TTT,nsimul))
loadingsIPCA_wrong=array(NA, dim = c(N,nfac0,TTT,nsimul))
gamma=array(NA, dim = c(nfacchar,nfac0,nsimul))
gammatilde=array(NA, dim = c(nfacchar,nfac0,nsimul))
gammatildewrong=array(NA, dim = c(1,nfac0,nsimul))
Htilde=array(NA, dim = c(nfac0,nfac0,TTT,nsimul))
common=array(NA,dim=c(TTT,N,nsimul))
commontildeIPCA=array(NA,dim=c(TTT,N,nsimul))
commontilde=array(NA,dim=c(TTT,N,nsimul))
loadingsIPCA2=array(NA, dim = c(N,nfac0,TTT,nsimul))
commontildeIPCAwrong=array(NA,dim=c(TTT,N,nsimul))

pe2=matrix(NA,TTT,nsimul)
pei2=matrix(NA,TTT,nsimul)
peiw2=matrix(NA,TTT,nsimul)
peo2=matrix(NA,TTT,nsimul)
peio2=matrix(NA,TTT,nsimul)
peiwo2=matrix(NA,TTT,nsimul)
r2=matrix(NA,1,nsimul)
r2i=matrix(NA,1,nsimul)
r2iw=matrix(NA,1,nsimul)
r2p=matrix(NA,1,nsimul)
r2ip=matrix(NA,1,nsimul)
r2iwp=matrix(NA,1,nsimul)
r2o=matrix(NA,1,nsimul)
r2io=matrix(NA,1,nsimul)
r2iwo=matrix(NA,1,nsimul)
data=array(NA,dim=c(TTT,N,nsimul))
HtildeIPCA=array(NA,dim=c(nfac0,nfac0,nsimul))


FtildeIPCAo=array(NA,dim=c(TTT,nfac0,nsimul))
FtildeIPCAowrong=array(NA,dim=c(TTT,nfac0,nsimul))
loadingstildeIPCAo=array(NA,dim=c(N,nfac0,TTT,nsimul))
loadingstildeIPCAowrong=array(NA,dim=c(N,nfac0,TTT,nsimul))
commontildeIPCAo=array(NA,dim=c(TTT,N,nsimul))
commontildeIPCAowrong=array(NA,dim=c(TTT,N,nsimul))

r2po=matrix(NA,1,nsimul)
r2ipo=matrix(NA,1,nsimul)
r2iwpo=matrix(NA,1,nsimul)



}



t_time = Sys.time()
supp = parLapply(cl = cl, X = 1:nsimul, f = Montecarlo4, nfacchar = nfacchar , factors0 = F0, seed0=seed0 ,
           q=q,   eps1=eps1,   indIPCA2=indIPCA2, indIPCA1=indIPCA1 ,    TT=TT  , hNT=hNT ,N = N, w = c("factors", "riskpremia", "sdf"))
Sys.time() - t_time
for (j in 1:nsimul)
{
 loadingsIPCA[,,,j]= supp[[j]][[1]]
  FtildeIPCA[,,j]=supp[[j]][[2]]
  FtildeIPCAwrong[,,j]=supp[[j]][[3]]
  loadingstildeIPCA[,,,j]=supp[[j]][[4]]
  loadingstildeIPCA_wrong[,,,j]=supp[[j]][[5]]
  loadingstilde[,,,j]=supp[[j]][[6]]
  Ftilde[,,j]=supp[[j]][[7]]
  gamma[,,j]=supp[[j]][[8]]
  gammatilde[,,j]=supp[[j]][[9]]
  gammatildewrong[,,j]=supp[[j]][[10]]
  Htilde[,,,j]=supp[[j]][[11]]
  common[,,j]=supp[[j]][[12]]#array(NA,dim=c(TTT,N,nsimul))
  commontildeIPCA[,,j]=supp[[j]][[13]]#array(NA,dim=c(TTT,N,nsimul))
  commontilde[,,j]=supp[[j]][[14]]#array(NA,dim=c(TTT,N,nsimul))
  loadingsIPCA2[,,,j]=supp[[j]][[15]]#array(NA, dim = c(N,nfac0,TTT,nsimul))
  commontildeIPCAwrong[,,j]=supp[[j]][[16]]#array(NA,dim=c(TTT,N,nsimul))
  pe2[,j]=supp[[j]][[17]]
  pei2[,j]=supp[[j]][[18]]
  peiw2[,j]=supp[[j]][[19]]
  r2[j]=supp[[j]][[20]]
  r2i[j]=supp[[j]][[21]]
  r2iw[j]=supp[[j]][[22]]
  r2p[j]=supp[[j]][[23]]
  r2ip[j]=supp[[j]][[24]]
  r2iwp[j]=supp[[j]][[25]]
  r2o[j]=supp[[j]][[26]]
  r2io[j]=supp[[j]][[27]]
  r2iwo[j]=supp[[j]][[28]]
  data[,,j]=supp[[j]][[29]]
  HtildeIPCA[,,j]=supp[[j]][[30]]
  
  FtildeIPCAo[,,j]=supp[[j]][[31]]
  FtildeIPCAowrong[,,j]=supp[[j]][[32]]
  loadingstildeIPCAo[,,,j]=supp[[j]][[33]]
  loadingstildeIPCAowrong[,,,j]=supp[[j]][[34]]
  commontildeIPCAo[,,j]=supp[[j]][[35]]
  commontildeIPCAowrong[,,j]=supp[[j]][[36]]
  
  r2po[j]=supp[[j]][[37]]
  r2ipo[j]=supp[[j]][[38]]
  r2iwpo[j]=supp[[j]][[39]]
  
  peo2[,j]=supp[[j]][[40]]
  peio2[,j]=supp[[j]][[41]]
  peiwo2[,j]=supp[[j]][[42]]
  
  
  
}

  

clt1=matrix(0,nsimul,1)
clt2=matrix(0,nsimul,1)
meanload0=matrix(0,nsimul,nfac0)
loadings0rot=array(NA,dim=c(N,nfac0,nsimul))
cr=0;
crl=0;
cri=0;
criw=0;
crli=0;
crliw=0;
crli2=0;
crliw2=0;
Adata=matrix(0,TTT,nsimul)
Ac=matrix(0,TTT,nsimul)
Actilde=matrix(0,TTT,nsimul)
ActildeIPCA=matrix(0,TTT,nsimul)


loadingstildeIPCA0=apply(loadingstildeIPCA,c(1,2,4),mean)
loadingstildeIPCA0_wrong=apply(loadingstildeIPCA_wrong,c(1,2,4),mean)
factors0rot=array(NA,dim=c(TTT,nfac0,nsimul))
loadings0rot=array(NA,dim=c(N,nfac0,TTT,nsimul))


crltime=0
crlitime=0
crliwtime=0
crlitime2=0
crliwtime2=0

loadingsIPCArot=array(NA,dim=c(N,nfac0,TTT,nsimul))
loadingsIPCA2rot=array(NA,dim=c(N,nfac0,TTT,nsimul))


for (j in 1:nsimul)
{
  
  
  t=TT
  for (t in 1:(TT-1))
  { loadingsIPCArot[,,t,j] =matrix(loadingsIPCA[,,t,j],N,nfac0)%*%solve(t(matrix(HtildeIPCA[,,j],nfac0,nfac0)))
  loadingsIPCA2rot[,,t,j] =matrix(loadingsIPCA2[,,t,j],N,nfac0)%*%solve(t(matrix(HtildeIPCA[,,j],nfac0,nfac0)))
  }
  
  
roll=TT
  for (roll in TT:TTT)
  {
    
    
   
      loadingsIPCArot[,,roll,j] =matrix(loadingsIPCA[,,roll,j],N,nfac0)%*%solve(t(matrix(HtildeIPCA[,,j],nfac0,nfac0)))
      loadingsIPCA2rot[,,roll,j] =matrix(loadingsIPCA2[,,roll,j],N,nfac0)%*%solve(t(matrix(HtildeIPCA[,,j],nfac0,nfac0)))
      
      
    
    
  F0rot=matrix(TT,nfac0)  
  F0rot=matrix(F0[(roll-TT+1):roll,],TT,nfac0)%*%Htilde[,,roll,j]
  factors0rot[roll,,j]=F0rot[TT,]
  #load0rot=matrix(0,N,nfac0)
 # load0rot=loadingstildeIPCA[,,(roll-TT+1):roll,j]%*%t(solve(Htilde[,,roll,j]))
 # loadings0rot[,,roll,j]=load0rot
  loadings0rot[,,roll,j]=matrix(apply(array(loadingsIPCA[,,(roll-TT+1):roll,j],dim=c(N,nfac0,TT,1)),c(1,2,4),mean),N,nfac0)%*%t(solve(Htilde[,,roll,j]))
  
  crltime=crltime+abs(cor(loadings0rot[,1,roll,j],loadingstilde[,1,roll,j]))
  crlitime=crlitime+abs(cor(loadingsIPCArot[,1,roll,j],loadingstildeIPCA[,1,roll,j]))
  crliwtime=crliwtime+abs(cor(loadingsIPCArot[,1,roll,j],loadingstildeIPCA_wrong[,1,roll,j]))
  crlitime2=crlitime2+abs(cor(loadingsIPCA2rot[,1,roll,j],loadingstildeIPCA[,1,roll,j]))
  crliwtime2=crliwtime2+abs(cor(loadingsIPCA2rot[,1,roll,j],loadingstildeIPCA_wrong[,1,roll,j]))
  
  
  
  }
  
F0rot=matrix(TT,nfac0)  
F0rot=matrix(F0[1:TT,],TT,nfac0)%*%Htilde[,,TT,j]
factors0rot[1:TT,,j]=F0rot



F0IPCArot=matrix(TTT,nfac0)  
F0IPCArot=matrix(F0,TTT,nfac0)%*%matrix(HtildeIPCA[,,j],nfac0,nfac0)




  cr=cr+abs(cor(factors0rot[,1,j],Ftilde[,1,j]))
  cri=cri+abs(cor(F0IPCArot[,1],FtildeIPCA[,1,j]))
  criw=criw+abs(cor(F0IPCArot[,1],FtildeIPCAwrong[,1,j]))
  
  
  # crltime=crltime+abs(cor(loadings0rot[,1,roll,j],loadingstilde[,1,roll,j]))
  # crlitime=crlitime+abs(cor(loadingsIPCArot[,1,roll,j],loadingstildeIPCA[,1,roll,j]))
  # crliwtime=crliwtime+abs(cor(loadingsIPCArot[,1,roll,j],loadingstildeIPCA_wrong[,1,roll,j]))
  # crlitime2=crlitime2+abs(cor(loadingsIPCA2rot[,1,roll,j],loadingstildeIPCA[,1,roll,j]))
  # crliwtime2=crliwtime2+abs(cor(loadingsIPCA2rot[,1,roll,j],loadingstildeIPCA_wrong[,1,roll,j]))
  
  # this is for the first time period only and first factor j=1?
  crl=crl+abs(cor(loadings0rot[,1,TT,j],loadingstilde[,1,TT,j]))
  crli=crli+abs(cor(loadingsIPCArot[,1,1,j],loadingstildeIPCA[,1,1,j]))
  
  crliw=crliw+abs(cor(loadingsIPCArot[,1,1,j],loadingstildeIPCA_wrong[,1,1,j]))
  
  crli2=crli2+abs(cor(loadingsIPCA2rot[,1,1,j],loadingstildeIPCA[,1,1,j]))
  
  crliw2=crliw2+abs(cor(loadingsIPCA2rot[,1,1,j],loadingstildeIPCA_wrong[,1,1,j]))
  
  Adata[1:(TTT-Tind[ii]+1),j]=eigen(data[Tind[ii]:TTT,,j]%*%t(data[Tind[ii]:TTT,,j]))$values
  Ac[1:(TTT-Tind[ii]+1),j]=eigen(common[Tind[ii]:TTT,,j]%*%t(common[Tind[ii]:TTT,,j]))$values
  Actilde[1:(TTT-Tind[ii]+1),j]=eigen(commontilde[Tind[ii]:TTT,,j]%*%t(commontilde[Tind[ii]:TTT,,j]))$values
  ActildeIPCA[1:(TTT-Tind[ii]+1),j]=eigen(commontildeIPCA[Tind[ii]:TTT,,j]%*%t(commontildeIPCA[Tind[ii]:TTT,,j]))$values
   
  
  #clt1[j]=sqrt(N/as.numeric(Ftilde_acm[1,1,1,j]) )*(Htilde[1,1,j]*F0[1,1]-Ftilde[1,1,j])
  #clt2[j]=sqrt(N/as.numeric(acmtildebai[1,1,1,j]) )*(Htilde[1,1,j]*F0[1,1]-Ftilde[1,1,j])
  #clt1[j]=sqrt(N/as.numeric(Ftilde_acm[1,1,1,j]) )*(F0rot[1,1]-Ftilde[1,1,j])
  #clt2[j]=sqrt(N/as.numeric(acmtildebai[1,1,1,j]) )*(F0rot[1,1]-Ftilde[1,1,j])
  
  
  #meanload0[j,]=as.numeric(apply(matrix(loadings0rot[,,j],N,nfac0),2,mean))
  #loadmeanIPCA[j]=
  #loadmeanIPCAwrong[j]
  
  
  
}
#loadmean1=t(as.matrix(loadmean[1,,],nfac0,N))

#cor(meanload0,t(loadmean1))
#cor(meanload0[,1],loadmean1[,1])


cr=cr/nsimul
#cr

crl=crl/nsimul  # this is time=TT
#crl

cri=cri/nsimul
#cri

crli=crli/nsimul
#crli

crli2=crli2/nsimul
#crli2


criw=criw/nsimul
#criw

crliw=crliw/nsimul
#crliw


crliw2=crliw2/nsimul
#crliw2


crlitime=crlitime/((TTT-TT+1)*nsimul)
#crlitime

crlitime2=crlitime2/((TTT-TT+1)*nsimul)
#crlitime2


crltime=crltime/((TTT-TT+1)*nsimul)
#crltime

crliwtime=crliwtime/((TTT-TT+1)*nsimul)
#crliwtime

crliwtime2=crliwtime2/((TTT-TT+1)*nsimul)
#crliwtime2



Adatamean=apply(Adata,1,mean)
Acmean=apply(Ac,1,mean)
Actildemean=apply(Actilde,1,mean)
ActildeIPCAmean=apply(ActildeIPCA,1,mean)

#cbind(mean(r2),mean(r2i),mean(r2iw))
#cbind(mean(r2p),mean(r2ip),mean(r2iwp))
#cbind(mean(r2o),mean(r2io),mean(r2iwo))
#cbind(mean(r2po),mean(r2ipo),mean(r2iwpo))


#cbind(cr,crl,cri,crli,crli2,criw,crliw,crliw2)
#cbind(crltime,crlitime,crlitime2,crliwtime,crliwtime2)


#cbind(Adatamean[1:10],Acmean[1:10],Actildemean[1:10],ActildeIPCAmean[1:10])

# jpeg(filename = "myplot.jpeg")
#  plot((1:1:TTT),log(Acmean),type="l")
#  lines(log(Actildemean),col=2)
#  lines(log(ActildeIPCAmean),col=3)
#  dev.off()

 jpeg(filename = paste("myplot_is_ipca1000_check","N",Nind[jj],"T",Tind[ii],"load",pca2[bb] , ".jpeg", sep=""))
  
 mm=max(apply(peiw2,1,mean))
   plot((1:(TTT-2*Tind[ii]+1)),apply(pe2[(2*Tind[ii]):TTT,],1,mean),ylim=c(0,mm),
        xlab="time",ylab= expression("pricing errors metric" ~ delta) , main=
          paste("Pricing Performance: In-Sample \nN=",Nind[jj],", T=",Tind[ii]))
#        "Pricing Performance: In-Sample")
   points((1:(TTT-2*Tind[ii]+1)),apply(pei2[(2*Tind[ii]):TTT,],1,mean),col=2)
   points((1:(TTT-2*Tind[ii]+1)),apply(peiw2[(2*Tind[ii]):TTT,],1,mean),col=3)
   legend("topleft", legend = c("local PCA", "IPCA", "IPCA (wrong)"), col = 1:3, lty = 1)      
 dev.off()
 

 jpeg(filename = paste("myplot_oos_ipca1000_check","N",Nind[jj],"T",Tind[ii],"load",pca2[bb] , ".jpeg", sep=""))
 mm=max(apply(peiwo2,1,mean))
 plot((1:(TTT-2*Tind[ii]+1)),apply(peo2[(2*Tind[ii]):TTT,],1,mean),ylim=c(0,mm),
      xlab="time",ylab= expression("pricing errors metric" ~ delta) , main=
        paste("Pricing Performance: Out-of-Sample \nN=",Nind[jj],", T=",Tind[ii]))
        #              "Pricing Performance: Out-of-Sample")
 #     xlab="time",ylab="pricing error", main="Pricing Performance")
 points((1:(TTT-2*Tind[ii]+1)),apply(peio2[(2*Tind[ii]):TTT,],1,mean),col=2)
 points((1:(TTT-2*Tind[ii]+1)),apply(peiwo2[(2*Tind[ii]):TTT,],1,mean),col=3)
 legend("topleft", legend = c("local PCA", "IPCA", "IPCA (wrong)"), col = 1:3, lty = 1)      
  dev.off()

 
# tables

# pricing errors

sink(file = "output_ipca1000check_2024.txt", append=TRUE)


print(cbind(Nind[jj],Tind[ii],indIPCA1,indIPCA2)) # when 1 loadings with rw; when 0 only IPCA part

print(cbind(mean(apply(pe2[(2*Tind[ii]):TTT,],1,mean)),mean(apply(pei2[(2*Tind[ii]):TTT,],1,mean)),
      mean(apply(peiw2[(2*Tind[ii]):TTT,],1,mean))))

print(cbind(mean(apply(peo2[(2*Tind[ii]):TTT,],1,mean)),mean(apply(peio2[(2*Tind[ii]):TTT,],1,mean)),mean(apply(peiwo2[(2*Tind[ii]):TTT,],1,mean))))


# maximum eigenvalue

print(cbind(Adatamean[1:10],Acmean[1:10],Actildemean[1:10],ActildeIPCAmean[1:10]))

# rsquare

print(cbind(mean(r2),mean(r2i),mean(r2iw)))
#cbind(mean(r2p),mean(r2ip),mean(r2iwp))
print(cbind(mean(r2o),mean(r2io),mean(r2iwo)))
#cbind(mean(r2po),mean(r2ipo),mean(r2iwpo))

# correlations

print(cbind(cr,crl,cri,crli,criw,crliw))
print(cbind(crltime,crlitime,crliwtime))

sink()

}
}
}
rm(supp, temp, j, t, t_time)

# T=10
# N=11
# F=matrix(rnorm(T^2,0,1),T,T)%*%rnorm(T,0,1)
# Lam=rnorm(N,0,1)
# X=F%*%t(Lam)+mcycleatrix(rnorm(T*N,0,1),T,N)
# 
# S=(X)%*%t(X)/(N*T)
# A=eigen(S)
# hatF=A$vectors[,1]*sqrt(T)
# hatLm=t(X)%*%hatF/T
# e=X - hatF%*%t(hatLm)
# t(e)%*%(hatF-matrix(1,T,1)*mean(hatF))/T
# mean((hatF-matrix(1,T,1)*mean(hatF)))









# here stop new Monte for revision - Oct 22.

# here old Monte first draft

# Factors : Compare ACM for a given t
t = TT
{
  #1: Population
  Ftilde_acm0[,,t]
  
  #2: mean of estimates across simulations
  method2_F = diag(0, nfac0)
  for (j in 1:nsimul) method2_F = method2_F + Ftilde_acm[,,t,j]/nsimul
  Ftilde_acm0[,,t]; method2_F
  
  #3: montecarlo 
  method3_F = matrix(Ftilde[t,,], ncol = nsimul)
  for (j in 1:nsimul) method3_F[,j] = sqrt(N)*(method3_F[,j] - t(Htilde[,,j]) %*% F0[t,])
  method3_F = method3_F %*% t(method3_F)/nsimul
  Ftilde_acm0[,,t]; method3_F
  
  rm(method2_F, method3_F, j)
}
# Compare ACM for every t = 1, ..., TT (ONLY FOR r = 1)
{
  if (nfac0 != 1) print ("This is only for r = 1")
  
  result = matrix(NA, ncol = 3, nrow = TT)
  for (t in 1:TT)
  {
    temp1 = Ftilde_acm0[,,t]
    
    supp = diag(0, nfac0)
    for (j in 1:nsimul) supp = supp + Ftilde_acm[,,t,j]
    temp2 = supp/nsimul
    
    temp = matrix(Ftilde[t,,], ncol = nsimul)
    for (j in 1:nsimul) temp[,j] = sqrt(N)*(temp[,j] - t(Htilde[,,j]) %*% F0[t,])
    temp3 = temp %*% t(temp)/nsimul
    
    result[t,] = c(temp1, temp2, temp3)
  }
  rm(temp, temp1, temp2, temp3, supp, t, j)
  
  result
  plot(x = 1:TT, y = result[,1], type = "l", ylim = c(min(result), max(result)))
  lines(x = 1:TT, y = result[,2], col = 2)
  lines(x = 1:TT, y = result[,3], col = 3)
}

# Risk premia
{
  #1: Population
  RiskPremiatilde_acm0
  
  #2: mean of estimates across simulations
  method2_RP = diag(0, nfac0)
  for (j in 1:nsimul) method2_RP = method2_RP + RiskPremiatilde_acm[,,j]/nsimul
  RiskPremiatilde_acm0; method2_RP
  
  #3: montecarlo
  method3_RP = matrix(RiskPremiatilde, ncol = nsimul)
  for (j in 1:nsimul) method3_RP[,j] = sqrt(N)*(method3_RP[,j] - t(Htilde[,,j]) %*% RiskPremiaPost)
  method3_RP = method3_RP %*% t(method3_RP)/nsimul
  RiskPremiatilde_acm0; method3_RP
  
  rm(method2_RP, method3_RP, j)
}
# SDF : Compare ACM for every t = 1, ..., TT
{
  acm_SDF = matrix(NA, ncol = 3, nrow = TT)
  acm_SDF[,1] = SDFtilde_acm0
  acm_SDF[,2] = apply(SDFtilde_acm, 1, mean)
  temp = sqrt(N) * (SDFtilde - SDFpost %*% t(rep(1, nsimul)))
  acm_SDF[,3] = apply(temp^2, 1, mean)
  
  acm_SDF
  plot(x = 1:TT, y = acm_SDF[,1], type = "l", ylim = c(min(0), max(acm_SDF)))
  lines(x = 1:TT, y = acm_SDF[,2], col = 2)
  lines(x = 1:TT, y = acm_SDF[,3], col = 3)
  
  rm(temp)
}



#other
{
plot(method3_F[4,])

t = TT
supp = matrix(Ftilde[t,,], ncol = nsimul)
for (j in 1:nsimul) supp[,j] = t(Htilde[,,j]) %*% F0[t,]

a = 1#4
plot(supp[a,])
lines(Ftilde[t, a,], type = "b", col = 2)
lines(rep((F0 %*% H)[t, a], nsimul), col = 3)
plot(Ftilde[t, a,]/supp[a,])

supp1 = matrix(Ftilde[t,,], ncol = nsimul)
for (j in 1:nsimul) supp1[,j] = sqrt(N)*(supp1[,j] - t(Htilde[,,j]) %*% F0[t,])
plot(supp1[a,])
mean(supp1[a,])

plot(sign(supp[a,]))
which(sign(supp[a,]) == -1)
which(sign(Ftilde[t, a,]) == -1)

t(F0) %*% F0 / TT ; det(t(F0) %*% F0 / TT)
t(Ftilde[,,1]) %*% Ftilde[,,1] / TT ; det(t(Ftilde[,,1]) %*% Ftilde[,,1] / TT)
t(F0 %*% Htilde[,,1]) %*% (F0 %*% Htilde[,,1]) / TT; det(t(F0 %*% Htilde[,,1]) %*% (F0 %*% Htilde[,,1]) / TT)
det(Htilde[,,1])^2 * det(t(F0) %*% F0 / TT)

t(Ftilde[,,1]) %*% F0 / TT; Qtilde[,,1]
t(Qtilde[,,1]) %*% sqrtm(Utilde_inv[,,1]); Psitilde[,,1]

Htilde[,,1]; det(Htilde[,,1])

supp = Htilde[,,1]
supp = Htilde[,,1] %*% sqrtm(Utilde_inv[,,1])
t(F0 %*% supp) %*% (F0 %*% supp) / TT; det(t(F0 %*% supp) %*% (F0 %*% supp) / TT)


supp = array(NA, dim = c(TT, nfac0, nsimul))
for (j in 1:nsimul) supp[,,j] = F0 %*% Htilde[,,j]
supp1 = array(NA, dim = c(nfac0, nfac0, nsimul))
for (j in 1:nsimul) supp1[,,j] = t(supp[,,j]) %*% supp[,,j]/TT
plot(supp1[1,1,]); mean(supp1[1,1,]); var(supp1[1,1,])
plot(supp1[2,2,]); mean(supp1[2,2,]); var(supp1[2,2,])

}


# CHECKS
{
# controllo 1: E coincida con B
{
  Ftilde_acm0_fun_check = function(t, factors0)
  {
    r = ncol(factors0)
    TT = nrow(factors0)
    
    Imatr = diag(1, r)
    ImatT = diag(1, TT)
    K_rT = commutation_matrix(r, TT)
    K_Tr = commutation_matrix(TT, r)
    
    temp = population_quantities(factors0)
    U_inv = temp[[1]]
    Psi = temp[[2]]
    Q = temp[[3]]
    H = temp[[4]]
    sigma2 = temp[[5]]
    sigma4 = temp[[6]]
    Sigma_lambda = temp[[7]]
    
    Ue = diag(sigma4, TT^2)
    for (j in 1:TT)
      for (k in 1:TT)
        Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
    
    E1 = (ImatT %x% (U_inv %*% (t(factors0 %*% H)/TT))) %*% Ue %*% (ImatT %x% ((factors0 %*% H/TT) %*% U_inv))
    E2 = sigma2*ImatT %x% U_inv
    E3 = (factors0 %*% Sigma_lambda %*% t(factors0)) %x% (sigma2/TT*U_inv^2)
    E4 = ((sigma2/TT* factors0 %*% H %*% U_inv) %x% t(factors0 %*% H)) %*% K_rT
    E5 = K_Tr %*% ((sigma2/TT*U_inv %*% t(factors0 %*% H)) %x% (factors0 %*% H))
    E = E1 + E2 + E3 + E4 + E5
    
    D = matrix(0, nrow = r, ncol = TT*r)
    D[, ((t - 1)*r + 1):(t*r)] = Imatr
    
    out = (D %*% E %*% t(D))
    
    return(out)
  }
  
  supp = matrix(0, nfac0, nfac0)
  for (t in 1:TT) 
  {
    print(t)
    print(Ftilde_acm0[,,t])
    print(Ftilde_acm0_fun_check(t, F0))
    supp = supp + abs(Ftilde_acm0[,,t] - Ftilde_acm0_fun_check(t, F0))
  }
  print(supp)
  rm(Ftilde_acm0_fun_check, t, supp)
}
# controllo 2: Etilde coincide con Btilde
{
Ftilde_acm_fun_check = function(t, k, factors, eigenvalue, valuef, sigma4)
{
  TT = nrow(factors)

  rf = rep(1, TT)
  
  IvecT = matrix(1, ncol = 1, nrow = TT)/TT
  ImatT = diag(1, TT)
  Imatk = diag(1, k)
  M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
  K_rT = commutation_matrix(k, TT)
  K_Tr = commutation_matrix(TT, k)
  
  gamma = apply(factors, 2, mean)
  omega = diag(1, k) - gamma %*% t(gamma)
  omega_inv = solve(omega)
  sigma2 = TT/(TT - k)*valuef
  U = diag(eigenvalue, k) - sigma2/TT*Imatk
  U_inv = solve(U)
  sigma_lambda = U
  Ue = diag(sigma4, TT^2)
  for (j in 1:TT)
    for (i in 1:TT)
      Ue[(j - 1)*TT + i, (i - 1)*TT + j] = Ue[(j - 1)*TT + i, (i - 1)*TT + j] + sigma4
  
  E1 = (ImatT %x% (U_inv %*% (t(factors)/TT))) %*% Ue %*% (ImatT %x% ((factors/TT) %*% U_inv))
  E2 = sigma2*ImatT %x% U_inv
  E3 = (factors %*% sigma_lambda %*% t(factors)) %x% (sigma2/TT*U_inv^2)
  E4 = ((sigma2/TT*factors %*% U_inv) %x% t(factors)) %*% K_rT
  E5 = K_Tr %*% ((sigma2/TT*U_inv %*% t(factors)) %x% factors)
  E = E1 + E2 + E3 + E4 + E5
  
  D = matrix(0, nrow = k, ncol = TT*k)
  D[, ((t - 1)*k + 1):(t*k)] = Imatk
  
  out = (D %*% E %*% t(D))
  
  return(out)
}
supp_fun = function(factors0, N)
{
  nfac0 = ncol(factors0)
  TT = nrow(factors0)
  
  # generazione dati
  loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  e = matrix(rnorm(TT*N, 0, 1), nrow = TT, ncol = N)
  
  common = factors0 %*% t(loadings)
  X = common + sqrt(nfac0)*e
  
  #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
  nfac_hat = nfac0
  
  temp = nfac_optimization(data = X, kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  valuef = temp[[1]]
  eigenv = temp[[2]]
  loadmean = temp[[3]]
  loadvar = temp[[4]]
  factors = as.matrix(temp[[5]])
  sigma4 = temp[[6]]
  gamma2pass = temp[[7]]
  rm(temp)
  
  valuef = valuef[nfac_hat]
  eigenv = eigenv[1:nfac_hat]
  factors = factors[,1:nfac_hat]
  sigma4 = sigma4[nfac_hat]
  
  sigma2 = TT/(TT - nfac_hat)*valuef
  Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
  Utilde_inv = solve(Utilde)
  Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  
  temp = population_quantities(factors0)
  Psi = temp[[2]]
  
  Ltilde = rep(1, nfac0)
  for (i in 1:nfac0) if (sign(Psi[1,i]) != sign(Psitilde[1,i])) Ltilde[i] = -Ltilde[i]
  
  factors = as.matrix(factors %*% diag(Ltilde))
  
  out1 = array(NA, dim = c(nfac_hat, nfac_hat, TT))
  out2 = array(NA, dim = c(nfac_hat, nfac_hat, TT))
  
  for (t in 1:TT) 
  {
    out1[,,t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
    out2[,,t] = Ftilde_acm_fun_check(t, nfac_hat, factors, eigenv, valuef, sigma4)
  }
  
  out = list(out1, out2)
  return(out)
}

supp = supp_fun(F0, N)
supp1 = matrix(0, nfac0, nfac0)
for (t in 1:TT) 
{
  print(t)
  print(supp[[1]][,,t])
  print(supp[[2]][,,t])
  supp1 = supp1 + abs(supp[[1]][,,t] - supp[[2]][,,t] )
}
print(supp1)

rm(Ftilde_acm_fun_check, supp_fun, t, supp, supp1)

}
# Controllo 3: Convergenza di Htilde ad H
{
  # vedere che H_tilde e H vicine a meno del segno per colonna
  for (i in 1:nfac0)
    for (j in 1:nfac0)
    {
      plot(Htilde[i, j, ])
      lines(rep(H[i, j], nsimul), col = 2)
    }
  temp = matrix(0, nfac0, nfac0)
  for (i in 1:nsimul) temp = temp + Htilde[,,i]/nsimul
  temp; H
  
  rm(i, j, temp)
}
# Controllo 4: Convergenza di Qtilde ad Q
{
  for (i in 1:nfac0)
    for (j in 1:nfac0)
    {
      plot(Qtilde[i, j, ])
      lines(rep(Q[i, j], nsimul), col = 2)
    }
  temp = matrix(0, nfac0, nfac0)
  for (i in 1:nsimul) temp = temp + Qtilde[,,i]/nsimul
  temp; Q
  
  rm(i, j, temp)
}
# Controllo 5: Ftilde%*%Htilde^-1 converge ad F?
{
  # Htilde ruota lo spazio da quello true spannato da F a quello tilde spannato da Ftilde.
  # Quindi per ogni simulazione, FHtilde giace su uno spazio diverso.
  # Tuttavia Htilde^-1 ruota lo spazio al contrario, da tilde a true. 
  # Quindi FtildeHtilde^-1 dovrebbero giacere tutti sullo stesso spazio true.
  # La media dei FtildeHtilde^-1 dovrebbe effere vicina a F
  temp = matrix(0, TT, nfac0)
  for (i in 1:nsimul) temp = temp + Ftilde[,,i] %*% solve(Htilde[,,i])/nsimul
  temp; F0
  
  rm(temp, i)
  # controllare se Ftilde converge a F%*%Htilde e' piu' complicato perche' si hanno valori diversi per ogni simulazione
}
# Controllo 6: F%*%Htilde converge ad F%*%H? 
{
  temp = matrix(0, TT, nfac0)
  for (i in 1:nsimul) temp = temp + F0 %*% Htilde[,,i]/nsimul
  F0 %*% H; temp
  
  rm(i, temp)
}
# Controllo 7: Convergenza di D_tilde a D
{
  Dtilde_check = function(t, k, factors, eigenvalue, valuef, sigma4)
  {
    TT = nrow(factors)

    rf = rep(1, TT)
    
    factors = as.matrix(factors[,1:k])
    eigenvalue = eigenvalue[1:k]
    valuef = valuef[k]
    sigma4 = sigma4[k]
    
    IvecT = matrix(1, ncol = 1, nrow = TT)/TT
    ImatT = diag(1, TT)
    Imatk = diag(1, k)
    M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
    K_rT = commutation_matrix(k, TT)
    K_Tr = commutation_matrix(TT, k)
    
    gamma = apply(factors, 2, mean)
    omega = diag(1, k) - gamma %*% t(gamma)
    omega_inv = solve(omega)
    sigma2 = TT/(TT - k)*valuef
    U = diag(eigenvalue, k) - sigma2/TT*Imatk
    U_inv = solve(U)
    sigma_lambda = U
    Ue = diag(sigma4, TT^2)
    for (j in 1:TT)
      for (k in 1:TT)
        Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
    
    D1 = (IvecT %x% Imatk) %*% omega_inv %*% (t(factors) %*% (ImatT[, t] - IvecT))
    D2 = ((ImatT[, t] - IvecT) %x% Imatk) %*% omega_inv %*% (t(factors) %*% IvecT)
    D31 = (1/TT * M_1T %*% factors %*% omega_inv %*% t(factors)) %x% (omega_inv %*% t(factors))
    D32 = (IvecT %x% (ImatT[, t] - IvecT)) + ((ImatT[, t] - IvecT) %x% IvecT)
    D = 1/rf[t]*(- D1 - D2 + D31 %*% D32)
    
    return(D)
  }
  D0_check = function(t, factors0)
  {
    TT = nrow(factors0)
    r = ncol(factors0)
    
    rf = rep(1, TT)
    
    IvecT = matrix(1, ncol = 1, nrow = TT)/TT
    ImatT = diag(1, TT)
    Imatk = diag(1, r)
    M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
    K_rT = commutation_matrix(r, TT)
    K_Tr = commutation_matrix(TT, r)
    
    temp = population_quantities(factors0)
    U_inv = temp[[1]]
    Psi = temp[[2]]
    Q = temp[[3]]
    H = temp[[4]]
    sigma2 = temp[[5]]
    sigma4 = temp[[6]]
    Sigma_lambda = temp[[7]]
    
    omega = (t(factors0) %*% M_1T %*% factors0)/TT
    omega_inv = solve(omega)
    Ue = diag(sigma4, TT^2)
    for (j in 1:TT)
      for (k in 1:TT)
        Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
    
    D1 = (IvecT %x% Imatk) %*% omega_inv %*% (t(factors0) %*% (ImatT[, t] - IvecT))
    D2 = ((ImatT[, t] - IvecT) %x% Imatk) %*% omega_inv %*% (t(factors0) %*% IvecT)
    D31 = (1/TT * M_1T %*% factors0 %*% omega_inv %*% t(factors0)) %x% (omega_inv %*% t(factors0))
    D32 = (IvecT %x% (ImatT[, t] - IvecT)) + ((ImatT[, t] - IvecT) %x% IvecT)
    D = (ImatT %x% solve(H)) %*% (1/rf[t]*(- D1 - D2 + D31 %*% D32)) 
    
    return(D)
  }
  supp_fun = function(factors0, N)
  {
    nfac0 = ncol(factors0)
    TT = nrow(factors0)
    
    # generazione dati
    loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
    e = matrix(rnorm(TT*N, 0, 1), nrow = TT, ncol = N)
    
    common = factors0 %*% t(loadings)
    X = common + sqrt(nfac0)*e
    
    #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
    nfac_hat = nfac0
    
    temp = nfac_optimization(data = X, kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
    valuef = temp[[1]]
    eigenv = temp[[2]]
    loadmean = temp[[3]]
    loadvar = temp[[4]]
    factors = as.matrix(temp[[5]])
    sigma4 = temp[[6]]
    gamma2pass = temp[[7]]
    rm(temp)
    
    valuef = valuef[nfac_hat]
    eigenv = eigenv[1:nfac_hat]
    factors = factors[,1:nfac_hat]
    sigma4 = sigma4[nfac_hat]
    
    sigma2 = TT/(TT - nfac_hat)*valuef
    Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
    Utilde_inv = solve(Utilde)
    Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
    
    temp = population_quantities(factors0)
    Psi = temp[[2]]

    Ltilde = rep(1, nfac0)
    for (i in 1:nfac0) if (sign(Psi[1,i]) != sign(Psitilde[1,i])) Ltilde[i] = -Ltilde[i]
    
    factors = as.matrix(factors %*% diag(Ltilde))
    
    Dtilde = matrix(NA, nrow = TT*nfac0, ncol = TT)
    for (t in 1:TT) Dtilde[,t] = Dtilde_check(t, nfac_hat, factors, eigenv, valuef, sigma4)
    
    return(Dtilde)
  }
  
  D0 = matrix(NA, nrow = TT*nfac0, ncol = TT)
  for (t in 1:TT) D0[,t] = D0_check(t, F0)
  Dtilde = array(NA, dim = c(TT*nfac0, TT, nsimul))
  for (i in 1:nsimul) Dtilde[,,i] = supp_fun(F0, N)
  
  t = 4
  supp = apply(Dtilde[,t,], 1, mean)
  cbind(supp, D0[,t])
  mean(supp - D0[,t])
  
  rm(Dtilde_check, D0_check, supp_fun, D0, Dtilde, t, i, supp)
}
# controllo 8: D vicino a derivata numerica di SDF componente per componente                             
{
  f1 = function(factors)
  {
    TT = nrow(factors)
    vecT = matrix(1, nrow = TT, ncol = 1)
    out = 1/TT * t(vecT) %*% factors
    return(out)
  }
  f2 = function(factors)
  {
    TT = nrow(factors)
    vecT = matrix(1, nrow = TT, ncol = 1)
    M1T = diag(1, TT) - vecT %*% solve(t(vecT) %*% vecT) %*% t(vecT)
    out = solve(1/TT * t(factors) %*% M1T %*% factors)
    return(out)
  }
  f3 = function(t, factors)
  {
    TT = nrow(factors)
    vecT = matrix(1, nrow = TT, ncol = 1)
    i_t  = matrix(0, nrow = TT, ncol = 1)
    i_t[t] = 1
    out = t(factors) %*% (i_t - vecT/TT)
    return(out)
  }
  D0_check = function(t, factors0)
  {
    TT = nrow(factors0)
    r = ncol(factors0)
    
    rf = rep(1, TT)
    
    IvecT = matrix(1, ncol = 1, nrow = TT)/TT
    ImatT = diag(1, TT)
    Imatk = diag(1, r)
    M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
    K_rT = commutation_matrix(r, TT)
    K_Tr = commutation_matrix(TT, r)
    
    temp = population_quantities(factors0)
    U_inv = temp[[1]]
    Psi = temp[[2]]
    Q = temp[[3]]
    H = temp[[4]]
    sigma2 = temp[[5]]
    sigma4 = temp[[6]]
    Sigma_lambda = temp[[7]]
    
    omega = (t(factors0) %*% M_1T %*% factors0)/TT
    omega_inv = solve(omega)
    Ue = diag(sigma4, TT^2)
    for (j in 1:TT)
      for (k in 1:TT)
        Ue[(j - 1)*TT + k, (k - 1)*TT + j] = Ue[(j - 1)*TT + k, (k - 1)*TT + j] + sigma4
    
    D1 = (IvecT %x% Imatk) %*% omega_inv %*% (t(factors0) %*% (ImatT[, t] - IvecT))
    D2 = ((ImatT[, t] - IvecT) %x% Imatk) %*% omega_inv %*% (t(factors0) %*% IvecT)
    D31 = (1/TT * M_1T %*% factors0 %*% omega_inv %*% t(factors0)) %x% (omega_inv %*% t(factors0))
    D32 = (IvecT %x% (ImatT[, t] - IvecT)) + ((ImatT[, t] - IvecT) %x% IvecT)

    D = (1/rf[t]*(- D1 - D2 + D31 %*% D32)) 
    
    return(list(-D1, -D2, D31 %*% D32, D))
  }
  
  t = 2
  SDFPost_fun(F0)[t]
  1 - f1(F0) %*% f2(F0) %*% f3(t, F0)
  dx = 0.001
  
  # component 1
  Jac1 = rep(NA, nfac0 *TT)
  for (j in 1:TT)
    for (i in 1:nfac0)
    {
      supp = matrix(0, nrow = TT, ncol = nfac0)
      supp[j, i] =  supp[j, i] + dx
      der = (f1(F0 + supp) - f1(F0 - supp))/(2*dx)
      Jac1[(j - 1) * nfac0 + i] = - der %*% f2(F0) %*% f3(t, F0)
    }
  cbind(D0_check(t, F0)[[1]], Jac1)
  
  # component 2
  Jac2 = rep(NA, nfac0 *TT)
  for (j in 1:TT)
    for (i in 1:nfac0)
    {
      supp = matrix(0, nrow = TT, ncol = nfac0)
      supp[j, i] =  supp[j, i] + dx
      der = (f2(F0 + supp) - f2(F0 - supp))/(2*dx)
      Jac2[(j - 1) * nfac0 + i] = - f1(F0) %*% der %*% f3(t, F0)
    }
  cbind(D0_check(t, F0)[[3]], Jac2)
  
  # component 3
  Jac3 = rep(NA, nfac0 *TT)
  for (j in 1:TT)
    for (i in 1:nfac0)
    {
      supp = matrix(0, nrow = TT, ncol = nfac0)
      supp[j, i] =  supp[j, i] + dx
      der = (f3(t, F0 + supp) - f3(t, F0 - supp))/(2*dx)
      Jac3[(j - 1) * nfac0 + i] = - f1(F0) %*% f2(F0) %*% der
    }
  cbind(D0_check(t, F0)[[2]], Jac3)
  
  # Tutto
  Jac = rep(NA, nfac0 *TT)
  for (j in 1:TT)
    for (i in 1:nfac0)
    {
      supp = matrix(0, nrow = TT, ncol = nfac0)
      supp[j, i] =  supp[j, i] + dx
      der = (SDFPost_fun(F0 + supp)[t] - SDFPost_fun(F0 - supp)[t])/(2*dx)
      Jac[(j - 1) * nfac0 + i] = der
    }
  cbind(D0_check(t, F0)[[4]], Jac1 + Jac2 + Jac3, Jac)
  
  rm(f1, f2, f3, D0_check, t, dx, Jac, Jac1, Jac2, Jac3, j, i, supp, der)
}
# Controllo 9: le identita'
{
  # identita' 1: t(H) %*% Sigma_F_hat %*% H = I_r
  # identita' 2: U_inv %*% solve(H) %*% Sigma_lambda %*% t(solve(H)) %*% U_inv = U_inv
  # identita' 3: U_inv %*% t(H) %*% Sigma_F_hat %*% Sigma_lambda = t(H)
  # identita' 4: (F %*% Sigma_lambda %*% t(F)) %x% (sigma2 * U_inv^2) = (F %*% Sigma_lambda %*% t(F)) %x% (sigma2 * U_inv %*% t(H) %*% Sigma_F_hat %*% H %*% U_inv)
  # ATTENZIONE: Alcune identita' non sono vere per ogni H, ma solo per H specifica. Motivo per il quale controlli 1 e 2 non e' detto che vengano.
  check = function(factors0)
  {
    TT = nrow(factors0)
    r = ncol(factors0)
    
    IvecT = matrix(1, ncol = 1, nrow = TT)/TT
    ImatT = diag(1, TT)
    Imatk = diag(1, r)
    M_1T = ImatT - IvecT %*% solve(t(IvecT) %*% IvecT) %*% t(IvecT)
    K_rT = commutation_matrix(r, TT)
    K_Tr = commutation_matrix(TT, r)
    
    temp = population_quantities(factors0)
    U_inv = temp[[1]]
    Psi = temp[[2]]
    Q = temp[[3]]
    H = temp[[4]]
    sigma2 = temp[[5]]
    sigma4 = temp[[6]]
    Sigma_lambda = temp[[7]]
    Sigma_F_hat = temp[[8]]
    
    omega = (t(factors0) %*% M_1T %*% factors0)/TT
    omega_inv = solve(omega)
   
    id1 = list(round(t(H) %*% Sigma_F_hat %*% H, 10), Imatk)
    id2 = list(round(U_inv %*% solve(H) %*% Sigma_lambda %*% t(solve(H)) %*% U_inv, 10), U_inv)
    id3 = list(U_inv %*% t(H) %*% Sigma_F_hat %*% Sigma_lambda, t(H))
    id4 = list(
      round((factors0 %*% Sigma_lambda %*% t(factors0)) %x% (sigma2 * U_inv^2), 10),
      round((factors0 %*% Sigma_lambda %*% t(factors0)) %x% (sigma2 * U_inv %*% t(H) %*% Sigma_F_hat %*% H %*% U_inv), 10)
    )
    return(list(id1, id2, id3, id4))
  }
  a = check(F0)
  sum(a[[1]][[1]] - a[[1]][[2]])
  sum(a[[2]][[1]] - a[[2]][[2]])
  sum(a[[3]][[1]] - a[[3]][[2]])
  sum(a[[4]][[1]] - a[[4]][[2]])
  
  rm(check, a)
}
# Controllo 10: controllare se Q continua a flippare segno
{  
  TT_fixed = TRUE
  Nmin = 100
  Nmax = 1000
  set.seed(2)
  
  if (TT_fixed == TRUE)
  {
    TT_i = TT
    F0_temp = matrix(NA, nrow = TT, ncol = nfac0)
    for (t in 1:TT) F0_temp[t,] = mvrnorm(1, meanF0, SigmaF0)
    Psi_temp = population_quantities(F0_temp)[[2]]
  }
  if (TT_fixed == FALSE)
  {
    F0_temp = matrix(NA, nrow = Nmax, ncol = nfac0)
    for (t in 1:Nmax) F0_temp[t,] = mvrnorm(1, meanF0, SigmaF0)
    Psi_all = array(NA, dim = c(nfac0, nfac0, Nmax))
  }
  
  Ltilde = matrix(NA, nrow = nfac0, ncol = Nmax)
  loadings = matrix(rnorm(Nmax*nfac0, 0, 1), nrow = Nmax, ncol = nfac0)
  
  for (N_i in Nmin:Nmax)
  { 
    print(N_i)
    
    if (TT_fixed == FALSE)
    {
      TT_i = N_i
      Psi_temp = population_quantities(F0_temp[1:TT_i,])[[2]]
      Psi_all[,,N_i] = Psi_temp
    }
    
    e = matrix(rnorm(TT_i*N_i, 0, 1), nrow = TT_i, ncol = N_i)
    common = F0_temp[1:TT_i,] %*% t(loadings[1:N_i,])
    X = common + sqrt(nfac0)*e
    
    temp = nfac_optimization(data = X, kmax = nfac0, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
    valuef = temp[[1]]
    eigenv = temp[[2]]
    loadmean = temp[[3]]
    loadvar = temp[[4]]
    factors = as.matrix(temp[[5]])
    sigma4 = temp[[6]]
    gamma2pass = temp[[7]]
    rm(temp)
    
    valuef = valuef[nfac0]
    eigenv = eigenv[1:nfac0]
    factors = factors[,1:nfac0]
    sigma4 = sigma4[nfac0]
    
    sigma2 = TT_i/(TT_i - nfac0)*valuef
    Utilde = diag(eigenv, nfac0) - sigma2/TT_i*diag(1, nfac0)
    Utilde_inv = solve(Utilde)
    Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT_i * t(F0_temp[1:TT_i,]) %*% factors)) %*% sqrtm((1/N_i * t(loadings[1:N_i,]) %*% loadings[1:N_i,])))
    
    Ltilde_supp = rep(1, nfac0)
    for (i in 1:nfac0) if (sign(Psi_temp[1,i]) != sign(Psitilde[1,i])) Ltilde_supp[i] = -Ltilde_supp[i]
    
#    factors = as.matrix(factors %*% diag(Ltilde_supp, nfac0))
#    Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT_i * t(F0_temp[1:TT_i,]) %*% factors)) %*% sqrtm((1/N_i * t(loadings[1:N_i,]) %*% loadings[1:N_i,])))
#    Ltilde_supp = rep(1, nfac0)
#    for (i in 1:nfac0) if (sign(Psi_temp[1,i]) != sign(Psitilde[1,i])) Ltilde_supp[i] = -Ltilde_supp[i]
    
    Ltilde[,N_i] = Ltilde_supp
    
    
  }
  plot(Ltilde[1,])
  plot(Ltilde[2,])
  Ltilde[,(Nmax - 100):Nmax]
  
  rm(TT_fixed, Nmin, Nmax, F0_temp, TT_i, N_i, t, i)
  rm(Ltilde, Ltilde_supp, Psi_temp, Psi_all)
  rm(loadings, e, common, X)
  rm(valuef, eigenv, loadmean, loadvar, factors, sigma4, gamma2pass, sigma2, Utilde, Utilde_inv, Psitilde)
}

}

#   montecarlo 3
# Input
N = 50
TT = 50
seed0=123
hNT=1/log(log(N))
nfac0 = 1
meanF0 = rep(0, nfac0)  
SigmaF0 = diag((1:nfac0), nfac0)
nsimul = 20#120

# Application starts
set.seed(2)
F0 = matrix(NA, nrow = TT, ncol = nfac0)
for (t in 1:TT) F0[t,] = mvrnorm(1, meanF0, SigmaF0)
#factors0=F0

# F0 = F0 - rep(1, TT) %*% t(apply(F0, 2, mean))    # delete these line at the end
# F0 = F0 %*% solve(chol(cov(F0) *(TT - 1)/ TT))    # delete these line at the end
# F0 = F0 %*% diag(sqrt(1:nfac0))                   # delete these line at the end
# t(F0) %*% F0 / TT                                 # delete these line at the end    

temp = population_quantities(F0)
U_inv = temp[[1]]
Psi = temp[[2]]
Q = temp[[3]]
H = temp[[4]]
sigma2 = temp[[5]]
sigma4 = temp[[6]]
Sigma_lambda = temp[[7]]

Ftilde_acm0 = array(NA, dim = c(nfac0, nfac0, TT))
for (t in 1:TT) Ftilde_acm0[,,t] = Ftilde_acm0_fun(t = t, factors = F0)

RiskPremiaPost = RiskPremiaPost_fun(F0)
RiskPremiatilde_acm0 = RiskPremiatilde_acm0_fun(F0)

SDFpost = SDFPost_fun(F0)
SDFtilde_acm0 = rep(NA, TT)
for (t in 1:TT) SDFtilde_acm0[t] = SDFtilde_acm0_fun(t, F0)

Qtilde = array(NA, dim = c(nfac0, nfac0, nsimul))
Htilde = array(NA, dim = c(nfac0, nfac0, nsimul))
Psitilde = array(NA, dim = c(nfac0, nfac0, nsimul))
Utilde_inv = array(NA, dim = c(nfac0, nfac0, nsimul))
Ftilde = array(NA, dim = c(TT, nfac0, nsimul))
loadingstilde = array(NA, dim = c(N, nfac0, nsimul))
loadings0 = array(NA, dim = c(N, nfac0, nsimul))

Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT, nsimul))
RiskPremiatilde = matrix(NA, nrow = nfac0, ncol = nsimul)
RiskPremiatilde_acm = array(NA, dim = c(nfac0, nfac0, nsimul))
SDFtilde = matrix(NA, nrow = TT, ncol = nsimul)
SDFtilde_acm = matrix(NA, nrow = TT, ncol = nsimul)
acmtildebai= array(NA, dim = c(nfac0, nfac0, TT, nsimul))
loadmean=array(NA, dim = c(nfac0, nfac0, nsimul)) 
FtildeIPCA = array(NA, dim = c(TT,nfac0,nsimul))
FtildeIPCAwrong = array(NA, dim = c(TT,nfac0,nsimul))
loadingstildeIPCA= array(NA, dim = c(N,nfac0,TT,nsimul))
loadingstildeIPCA_wrong=array(NA, dim = c(N,nfac0,TT,nsimul))
loadingsIPCA= array(NA, dim = c(N,nfac0,TT,nsimul))
loadingsIPCA_wrong=array(NA, dim = c(N,nfac0,TT,nsimul))






#set.seed(123)
t_time = Sys.time()
supp = parLapply(cl = cl, X = 1:nsimul, f = Montecarlo3, factors0 = F0, seed0=seed0 
                 , hNT=hNT ,N = N, w = c("factors", "riskpremia", "sdf"))
Sys.time() - t_time
for (j in 1:nsimul)
{
  Qtilde[,,j] = supp[[j]][[1]]
  Htilde[,,j] = supp[[j]][[2]]
  Psitilde[,,j] = supp[[j]][[3]]
  Utilde_inv[,,j] = supp[[j]][[4]]
  
  Ftilde[,,j] = supp[[j]][[5]]
  Ftilde_acm[,,,j] = supp[[j]][[6]]
  
  RiskPremiatilde[,j] = supp[[j]][[7]]
  RiskPremiatilde_acm[,,j] = supp[[j]][[8]]
  
  SDFtilde[,j] = supp[[j]][[9]]
  SDFtilde_acm[,j] = supp[[j]][[10]]
  loadingstilde[,,j]=supp[[j]][[11]]
  loadings0[,,j]=supp[[j]][[12]]
  
  acmtildebai[,,,j]=supp[[j]][[13]]
  loadmean[,,j]=supp[[j]][[14]]
  
  # comment below when montecarlo and montecarlo2
  FtildeIPCA[,,j] = supp[[j]][[15]]
  
  FtildeIPCAwrong[,,j] = supp[[j]][[16]]
  loadingstildeIPCA[,,,j]= supp[[j]][[17]]
  loadingstildeIPCA_wrong[,,,j]=supp[[j]][[18]]
  
  loadingsIPCA[,,,j]= supp[[j]][[19]]
  loadingsIPCA_wrong[,,,j]=supp[[j]][[20]]
  
  
}
#out = list(Qtilde1, Htilde2, Psitilde3, Utilde_inv4, factors5, Ftilde_acm6, RiskPremiatilde7, RiskPremiatilde_acm8, SDFtilde9, SDFtilde_acm10,loadingshat11,
#           loadings12,Ftilde_acmbai13,loadmean14,FtildeIPCA15,FtildeIPCA_wrong16,
#           loadingstildeIPCA17,loadingstildeIPCA_wrong18,loadingsIPCA19,loadingsIPCA_wrong20)

F0[1,1]
Ftilde[1,1,1]
loadings0[1,1,1]

# corr F and tildeF

clt1=matrix(0,nsimul,1)
clt2=matrix(0,nsimul,1)
meanload0=matrix(0,nsimul,nfac0)
loadings0rot=array(NA,dim=c(N,nfac0,nsimul))
cr=0;
crl=0;
cri=0;
criw=0;
crli=0;
crliw=0;

loadingstildeIPCA0=apply(loadingstildeIPCA,c(1,2,4),mean)
loadingstildeIPCA0_wrong=apply(loadingstildeIPCA_wrong,c(1,2,4),mean)


for (j in 1:nsimul)
{
  #  multiply byHj then get entry TT
  
  F0rot=F0%*%Htilde[,,j]
  loadings0rot[,,j]=loadings0[,,j]%*%t(solve(Htilde[,,j]))
  cr=cr+abs(cor(F0rot[,1],Ftilde[,1,j]))
  crl=crl+abs(cor(loadings0rot[,1,j],loadingstilde[,1,j]))
  
  cri=cri+abs(cor(F0[,1],FtildeIPCA[,1,j]))
  criw=criw+abs(cor(F0[,1],FtildeIPCAwrong[,1,j]))
  crli=crli+abs(cor(loadings0[,1,j],loadingstildeIPCA0[,1,j]))
  crliw=crliw+abs(cor(loadings0[,1,j],loadingstildeIPCA0_wrong[,1,j]))
  
  
  
  #clt1[j]=sqrt(N/as.numeric(Ftilde_acm[1,1,1,j]) )*(Htilde[1,1,j]*F0[1,1]-Ftilde[1,1,j])
  #clt2[j]=sqrt(N/as.numeric(acmtildebai[1,1,1,j]) )*(Htilde[1,1,j]*F0[1,1]-Ftilde[1,1,j])
  clt1[j]=sqrt(N/as.numeric(Ftilde_acm[1,1,1,j]) )*(F0rot[1,1]-Ftilde[1,1,j])
  clt2[j]=sqrt(N/as.numeric(acmtildebai[1,1,1,j]) )*(F0rot[1,1]-Ftilde[1,1,j])
  
  
  meanload0[j,]=as.numeric(apply(matrix(loadings0rot[,,j],N,nfac0),2,mean))
  #loadmeanIPCA[j]=
  #loadmeanIPCAwrong[j]
  
}
loadmean1=t(as.matrix(loadmean[1,,],nfac0,N))

cor(meanload0,t(loadmean1))
cor(meanload0[,1],loadmean1[,1])


cr=cr/nsimul
cr

crl=crl/nsimul
crl

cri=cri/nsimul
cri

crli=crli/nsimul
crli

criw=criw/nsimul
criw

crliw=crliw/nsimul
crliw



#  CLT of factors estimates
hist(clt1)
hist(clt2)

clt1=as.matrix(clt1)
clt2=as.matrix(clt2)


clt1h <- hist(clt1 , breaks = 30 , freq=FALSE, plot = FALSE)
clt2h <- hist(clt2 , breaks = 30 , freq=FALSE, plot = FALSE)


#hgBt <- hist(t_gro[,1] , breaks = 30 , plot = FALSE) # Save 2nd histogram data


# construct new colors - alpha gives the strength of the color
c1 <- rgb(173,216,230, max = 255, alpha = 80, names = "lt.blue")
c2 <- rgb(255,192,203, max = 255, alpha = 30, names = "lt.pink")
c3 <- rgb(144,238,144, max = 255, alpha = 80, names = "lt.green")
# c4 <- rgb(1,1,255, max = 255, alpha = 80, names = "lt.bo")

#  
#pdf('canadah2.pdf')
plot(clt1h , freq=FALSE,main="HISTOGRAM and PDF: \nDREFUS PORTFOLIO RETURN" ,xlim=c(-6,6),
     xlab="Monthly Returns",col = c1) # Plot 1st histogram using a transparent color

curve(dnorm(x, mean=mean(clt1), sd=sd(clt1)), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")
# add another asset
plot(clt2h,  col = c3, freq=FALSE,add = TRUE) # Add 2nd histogram using different color

#  START: now over long samples of overlapping rolling windows
#  but for IPCA over full samples

# IPCA over full sample TTT + localPCA

# # just IPCA  
# Montecarlo4i = function(i, seed0 ,factors0, TT, hNT, N, w = c("factors", "riskpremia", "sdf"))
{
  #   
  #   # this part including estimation IPCA does not roll
  #   
  #   #set.seed(seed0)
  #   nfac0 = ncol(factors0) # long time series TTT
  #   TTT = nrow(factors0)
  #   TT=TT
  #   
  #   
  #   
  #   sigmaz=1
  #   loadings = matrix(rnorm(N*nfac0, 0, 1), nrow = N, ncol = nfac0)
  #   t=2
  #   xi=matrix(0,(N*nfac0),TTT)
  #   for (t in 2:TTT)
  #   {
  #     
  #     xi[,t]=xi[,t-1]+rnorm((N*nfac0), 0, sigmaz)
  #     
  #     
  #   }
  #   
  #   
  #   
  #   e = matrix(rnorm(TTT*N, 0, 1), nrow = TTT, ncol = N)
  #   
  #   #  loadings = loadings - rep(1, N) %*% t(apply(loadings, 2, mean))     # delete these line at the end
  #   #  loadings = loadings %*% solve(chol(cov(loadings) *(N - 1)/ N))      # delete these line at the end
  #   #  t(loadings) %*% loadings / N                                        # delete these line at the end
  #   
  #   #common = factors0 %*% t(loadings)
  #   
  #   hNT=hNT #1/log(log(N))
  #   #hNT=10
  #   
  #   # generazione dati
  #   #  number char = 2
  #   nfacchar=2
  #   TTT1=2*TTT
  #   t=2
  #   meanchar=rnorm(nfacchar,0,1)#zero?
  #   char=matrix(0,(N*nfacchar),TTT1)
  #   for (t in 2:TTT1) 
  #   {
  #     char[,t]=(1-0.8)*(meanchar%x%matrix(1,N,1))  +  0.8*char[,t-1]+rnorm(N*nfacchar,0,1)
  #   }  
  #   
  #   gammab1=matrix(rnorm((nfac0*nfacchar),0,1),nfacchar,nfac0)
  #   #gammab0=rnorm((nfac0*nfacchar),0,1)
  #   
  #   charfinal=char[,(TTT1-TTT+1):TTT1]
  #   
  #   # ortho char
  #   orthchar=array(NA, dim = c(nfacchar, N , TTT1))
  #   t=2
  #   for (t in 2:TTT1)
  #   {
  #     mchar=matrix(char[,t],nfacchar,N)
  #     orthchar[,,t]=solve(sqrtm(mchar%*%t(mchar)))%*%mchar
  #   }
  #   
  #   
  #   #   orthchar[,,t]%*%t(orthchar[,,t])
  #   
  #   common=matrix(0,TTT,N)
  #   loadingsIPCA=array(0,dim=c(N,nfac0,TTT))
  #   loadingsIPCA2=array(0,dim=c(N,nfac0,TTT))
  #   loadingsIPCA_wrong=array(0,dim=c(N,nfac0,TTT))
  #   
  #   
  #   # common[1,]=(factors0[1,])%*%
  #   #   (t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #   #      hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)])
  #   # 
  #   # loadingsIPCA[,,1]=as.matrix(t(t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #   #                                 hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)]),N,nfac0)
  #   
  #   common[1,]=(factors0[1,])%*%
  #     (
  #    #   t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #        hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)])
  #   
  #   loadingsIPCA[,,1]=as.matrix(t(t(loadings + hNT*matrix(xi[,1] , nrow = N, ncol = nfac0)) +
  #                                   hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT)]),N,nfac0)
  #    t=1
  #   loadingsIPCA2[,,t]=t( hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)]) 
  #   
  #   t=2
  #   for (t in 2:TTT)
  #   {
  #     
  #     
  #     # common[t,]=(factors0[t,])%*%
  #     #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
  #     #      hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)])
  #      
  #     common[t,]=(factors0[t,])%*%
  #       (
  #         #t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
  #          hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)])
  #     
  #       
  #     loadingsIPCA[,,t]=t(    (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
  #                                hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)]) )
  #     
  #     loadingsIPCA2[,,t]=t( hNT*t(gammab1)%*%orthchar[,,(TTT1-TTT+t-1)]) 
  #     
  #   }
  #   
  #   X = common + sqrt(nfac0)*e
  #   X=as.matrix(X,TTT,N)
  #   
  #   #  nfac_hat = nfac_optimization(data = X, kmax = kmax, eps = eps, std = TRUE)[1]
  #   
  #   #  estimates IPCA
  #   commontildeIPCA=matrix(0,TTT,N)
  #   USU=matrix(0,nfacchar,nfacchar)
  #   tt=1
  #   for (tt in 1:TTT)
  #   {
  #     USU=USU +   t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)]))%*%X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])
  #   }
  #   
  #   USU.svd=svd(USU)
  #   
  #   
  #   aaa= as.matrix(USU.svd$u,nfacchar,nfacchar)
  #   gammab1tilde=aaa[,1:nfac0] #  kc times k
  #   gammab1tilde_wrong=aaa[1,1:nfac0] #  kcwrong=1 times k
  #   
  #   
  #   FtildeIPCA=matrix(0,TTT,nfac0)
  #   FtildeIPCA_wrong=matrix(0,TTT,nfac0)
  #   
  #   tt=1
  #   for (tt in 1:TTT)
  #   {
  #     FtildeIPCA[tt,]=t(X[tt,]%*%t(orthchar[,,(TTT1-TTT+tt-1)])%*%gammab1tilde) 
  #     FtildeIPCA_wrong[tt,]=t(X[tt,]%*%t(matrix(orthchar[1,,(TTT1-TTT+tt-1)],1,N))%*%gammab1tilde_wrong) 
  #     
  #     
  #     # commontildeIPCA[tt,]=t(FtildeIPCA[tt,])%*%
  #     #   (t(loadings + hNT*matrix(xi[,t] , nrow = N, ncol = nfac0)) +
  #     #      hNT*t(gammab1)%*%orthchar[,,(TTT-TT+t-1)])
  #     # 
  #     
  #     #t(gammab1tilde)%*%orthchar[,,(TTT-TT+tt-1)]%*%t(X[tt,])
  #   }
  #   
  #   # # this is alternative estimator for gammab
  #   # 
  #   # st=1;
  #   # den=matrix(0,2*nfac0,2*nfac0)
  #   # num=matrix(0,2*nfac0,1)
  #   # for (t in 1:TTT-1)
  #   # {
  #   #   den=den+ kronecker((orthchar[,,(TTT1-TTT+t)]%*%t(orthchar[,,(TTT1-TTT+t)])),FtildeIPCA[t+1,]%*%t(FtildeIPCA[t+1,]) )
  #   #   num=num+ (kronecker(orthchar[,,(TTT1-TTT+t)],FtildeIPCA[t+1,]))%*%X[t+1,]
  #   # }
  #   # 
  #   # gammabtilde2=matrix(solve(den)%*%num,2,nfac0)
  #   
  #   loadingstildeIPCA=array(0,dim=c(N,nfac0,TTT))
  #   loadingstildeIPCA_wrong=array(0,dim=c(N,nfac0,TTT))
  #   commonstildeIPCA=array(0,dim=c(TTT,N))
  #   commonstildeIPCAwrong=array(0,dim=c(TTT,N))
  #   
  #   
  #   
  #   t=1
  #   for (t in 1:TTT)
  #   {
  #     loadingstildeIPCA[,,t]=t(orthchar[,,(TTT1-TTT+t-1)])%*%gammab1tilde 
  #     loadingstildeIPCA_wrong[,,t]=as.matrix(orthchar[1,,(TTT1-TTT+t-1)],N,1)*gammab1tilde_wrong
  #     commontildeIPCA[t,]=FtildeIPCA[t,]%*%t(as.matrix(loadingstildeIPCA[,,t],N,nfac0))
  #     commontildeIPCAwrong[t,]=FtildeIPCA[t,]%*%t(as.matrix(loadingstildeIPCA_wrong[,,t],N,nfac0))
  #     
  #       #t(gammab1tilde)%*%orthchar[,,(TTT-TT+tt-1)]%*%t(X[tt,])
  #   }
  #   
  #   # find HtildeIPCA? not need as gammab is fixed
  #   
  #   # rolling estimation of local PCA
  #  #  comment PCA from here  to test IPCA
  # 
  #   nfac_hat = nfac0  
  #   
  #   RO=TTT-TT+1
  #   #TT = nrow(X)
  #   #N = ncol(X)
  #   # eigenv = matrix(NA,nfac_hat,RO) #eigen(1/(TT*N) * X %*% t(X))$values[1:nk]
  #   # valuef = array(NA,dim=c(nfac_hat,RO)) #rep(NA, nk)
  #   #                m_loadings =matrix(0, ncol = nk, nrow = nk)
  #   #                s_loadings =matrix(NA, ncol = nk, nrow = nk)
  #   factors = array(NA,dim=c(TT,nfac_hat,TTT))#matrix(NA, ncol = nk, nrow = TT)
  #   
  #   # sigma_4 = matrix(NA,nfac_hat,RO)#rep(NA, nk)
  #   # gamma_twopass = array(NA,dim=c(nfac_hat,nfac_hat,RO)) #matrix(NA, ncol = nk, nrow = nk)
  #   loadingshat = array(NA,dim=c(N,nfac_hat,TTT))#matrix(0, ncol = nk, nrow = nk)
  #   factorshat = array(NA,dim=c(TTT,nfac_hat))#matrix(0, ncol = nk, nrow = nk)
  #   Htilde=array(NA,dim=c(nfac0,nfac0,TTT))
  #   commontilde = matrix(NA,TTT,N)#matrix(0, ncol = nk, nrow = nk)
  #   
  #   {
  #   roll=TT
  #   for (roll in TT:TTT)
  #   {
  # 
  #     #temp = nfac_optimization(data = X[(roll-TT+1):roll,], kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  # 
  #     #  TT = nrow(X)
  #     #  N = ncol(X)
  #     nk=nfac_hat
  #     eig = eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$values[1:nk]
  # 
  #     val = rep(NA, nk)
  #     m_loadings = matrix(0, ncol = nk, nrow = nk)
  #     s_loadings = matrix(NA, ncol = nk, nrow = nk)
  #     #factors = matrix(NA, ncol = nk, nrow = TT)
  #     sigma_4 = rep(NA, nk)
  #     gamma_twopass = matrix(NA, ncol = nk, nrow = nk)
  # 
  # 
  #     for (k in 1:nk)
  #     {
  #       temp = sqrt(TT) * as.matrix(eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$vectors[,1:k])
  #       loadings = 1/TT *(t(X[(roll-TT+1):roll,]) %*% temp)
  #       m_loadings[1:k, k] = apply(loadings, 2, mean)
  #       val[k] = 1/(N*TT) * sum((X[(roll-TT+1):roll,] - temp %*% t(loadings))^2)
  #       temp1 = sum(diag((temp %*% t(temp))/TT)^2)
  #       temp1 = 3 + 27/TT*temp1 + 18*k/TT
  #       temp2 = 1/(N*TT) * sum((X[(roll-TT+1):roll,] - temp %*% t(loadings))^4)
  #       sigma_4[k] = temp2/temp1
  #       gamma_twopass[1:k,k] = solve(t(loadings) %*% loadings) %*% t(loadings) %*% apply(X[(roll-TT+1):roll,], 2, mean)
  # 
  #       #factors[1:TT, 1:k] = temp
  #       s_loadings[1:k, 1:k] = (t(loadings) %*% loadings)/N
  #     }
  # 
  #     valuef = val[nfac_hat]
  #     eigenv = eig[1:nfac_hat]
  #     s_loadings=s_loadings[1:nfac_hat,1:nfac_hat]
  #     gamma_twopass=gamma_twopass[1:nfac_hat,1:nfac_hat]
  #     m_loadings=m_loadings[1:nfac_hat,1:nfac_hat]
  # 
  #     #temp = nfac_optimization(data = X[roll-TT+1:roll,], kmax = nfac_hat, eps = c(1e-8, 0.175), std = F, est_extract = TRUE)
  #     # valuef = temp[[1]]
  #     # eigenv = temp[[2]]
  #     # loadmean= temp[[3]]
  #     # loadvar = temp[[4]]
  #     #factors = as.matrix(temp[[5]])
  #     factors[,,roll] =  sqrt(TT) * as.matrix(eigen(1/(TT*N) * X[(roll-TT+1):roll,] %*% t(X[(roll-TT+1):roll,]))$vectors[,1:nfac_hat])
  # 
  # 
  #     #sigma4 = temp[[6]]
  #     loadingshat[,,roll]=t(X[(roll-TT+1):roll,])%*%factors[,,roll]/TT
  #     #gamma2pass= temp[[7]]
  #     # rm(temp)
  # 
  #     #valuef = valuef[nfac_hat]
  #     #eigenv = eigenv[1:nfac_hat]
  #     # factors = as.matrix(factors[,1:nfac_hat])
  #     #  loadingshat = as.matrix(loadingshat[,1:nfac_hat])
  #     #  loadings = as.matrix(loadings[,1:nfac_hat])
  # 
  #     sigma_4 = sigma_4[nfac_hat]
  #     #
  #     sigma2 = TT/(TT - nfac_hat)*valuef
  #     Utilde = diag(eigenv, nfac_hat) - sigma2/TT*diag(1, nfac_hat)
  #     Utilde_inv = solve(Utilde)
  #     Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0[(roll-TT+1):roll,]) %*% factors[,,roll])) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  #     #
  #     temp = population_quantities(factors0[(roll-TT+1):roll,])
  #     Psi = temp[[2]]
  # 
  #     Ltilde = rep(1, nfac0)
  #     for (i in 1:nfac0)
  #     {
  #       t = 1
  #       while ((sign(round(Psi[t, i], 10)) == 0) | (sign(round(Psitilde[t,i], 10)) == 0)) t = t + 1
  #       if (sign(Psi[t,i]) != sign(Psitilde[t,i])) Ltilde[i] = -Ltilde[i]
  #     }
  # 
  #     factors[,,roll] = as.matrix(factors[,,roll] %*% diag(Ltilde, nfac0))
  #     loadingshat[,,roll] = as.matrix(loadingshat[,,roll] %*% diag(Ltilde, nfac0))
  # 
  #     factorshat[roll,]=factors[TT,,roll]
  # 
  #     meanloadings0=apply(loadingsIPCA[,,(roll-TT+1):roll],c(1,2),mean)
  #     loadings=meanloadings0
  #     # Qtilde[,,roll] = t((1/TT * t(factors0[(roll-TT+1):roll,]) %*% factors[,,roll]))
  #     Htilde[,,roll] = (1/N * t(loadings) %*% loadings)  %*% (1/TT * t(factors0[(roll-TT+1):roll,]) %*% factors[,,roll]) %*% Utilde_inv
  # 
  # 
  # 
  #     # Psitilde = t(sqrtm(Utilde_inv) %*% t((1/TT * t(factors0) %*% factors)) %*% sqrtm((1/N * t(loadings) %*% loadings)))
  #     #
  #     # # factors
  #     # Ftilde_acm = array(NA, dim = c(nfac0, nfac0, TT))
  #     # Ftilde_acmbai = array(NA, dim = c(nfac0, nfac0, TT))
  #     #
  #     # if (any(w == "factors"))
  #     # {
  #     #   for (t in 1:TT) Ftilde_acm[, , t] = Ftilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  #     #   for (t in 1:TT) Ftilde_acmbai[, , t] = Ftilde_acm_funbai(t, nfac_hat, factors, eigenv, valuef, sigma4)
  #     #
  #     # }
  #     # # risk premia
  #     # RiskPremiatilde = rep(NA, nfac0)
  #     # RiskPremiatilde_acm = matrix(NA, nfac0, nfac0)
  #     # if (any(w == "riskpremia"))
  #     # {
  #     #   RiskPremiatilde = RiskPremiatilde_fun(nfac_hat, factors)
  #     #   RiskPremiatilde_acm = RiskPremiatilde_acm_fun(nfac_hat, factors, eigenv, valuef, sigma4)
  #     # }
  #     # # SDF
  #     # SDFtilde = rep(NA, TT)
  #     # SDFtilde_acm = rep(NA, TT)
  #     # if (any(w == "sdf"))
  #     # {
  #     #   SDFtilde = SDFtilde_fun(nfac_hat, factors)
  #     #   for (t in 1:TT) SDFtilde_acm[t] = SDFtilde_acm_fun(t, nfac_hat, factors, eigenv, valuef, sigma4)
  #     # }
  #   }
  # 
  #   factorshat[1:TT,]=factors[,,TT]
  #   
  #   tt=1
  #   for (tt in 1:TTT)
  #   {
  #   commontilde[tt,]=t(factorshat[tt,])%*%t(loadingshat[,,tt])
  # }
  #   }
  #   # out = list(Qtilde, Htilde, Psitilde, Utilde_inv, factors, Ftilde_acm, RiskPremiatilde, RiskPremiatilde_acm, SDFtilde, SDFtilde_acm,loadingshat,
  #   #            loadings,Ftilde_acmbai,loadmean,FtildeIPCA,FtildeIPCA_wrong,
  #   #            loadingstildeIPCA,loadingstildeIPCA_wrong,loadingsIPCA,loadingsIPCA_wrong)
  #   out = list(loadingsIPCA,FtildeIPCA,FtildeIPCA_wrong,loadingstildeIPCA,loadingstildeIPCA_wrong,loadingshat,factorshat,gammab1,gammab1tilde,
  #              gammab1tilde_wrong,Htilde,common,commontildeIPCA,commontilde,loadingsIPCA2,commontildeIPCAwrong)
  #   
  #   return(out)
}


# canonical corr
n=100
X=matrix(rnorm(n*2,0,1),n,2)
Y=matrix(rnorm(n*2,0,1),n,2)
Z=matrix(rnorm(n*2,0,1),n,2)
W=matrix(rnorm(n*2,0,1),n,2)


#Y=cbind(X[,1:3])
cxy=cancor(X,Y)
cxy

 
# RGCCA

library(RGCCA)



#Canonical Correlation Analysis
cca.with.rgcca = rgcca(A= list(X, Y),
                       C = matrix(c(0, 1, 1, 0), 2, 2),
                       tau = c(0, 0))


# Generalized CCA

# special case J=2
J=2
C = matrix(c(0, 0, 1,
             0, 0, 1,
             1, 1,   0), J+1, J+1)
#aaa=rgcca(A= list(X,Y,Z))
gcca.with.rgcca = rgcca(A= list(X,Y,cbind(X,Y)), 
                        C = C,
                        tau = rep(0, J+1), scheme = "factorial")




# special case J=3
J=3
C = matrix(c(0, 0,  0,  1,
             0, 0,  0, 1,
             0, 0,  0, 1,
             1, 1,  1, 0), J+1, J+1)
#gcca.with.rgcca = rgcca(A= list(X,Y,Z,cbind(X,Y,Z)), 
gcca.with.rgcca = rgcca(A= list(X,Y,Z,W), 

                        C = C,
                        tau = rep(0, J+1), scheme = "factorial")
# J=4
J=4
C = matrix(c(0, 0,  0, 0, 1,
             0, 0,  0, 0, 1,
             0, 0,  0, 0, 1,
             0, 0,  0, 0, 1,
             1, 1,  1, 1, 0), J+1, J+1)
gcca.with.rgcca = rgcca(A= list(X,Y,Z,W,X), 
                        #gcca.with.rgcca = rgcca(A= list(X,Y,Z,W), 
                        
                        C = C,
                        tau = rep(0, J+1), scheme = "factorial")

# J=3 modified
J=3
C = matrix(c(0, 1,   1,
             1, 0, 1,
             1, 1,  0 ), J, J)
#gcca.with.rgcca = rgcca(A= list(X,Y,Z,cbind(X,Y,Z)), 
gcca.with.rgcca = rgcca(A= list(X,Y,Z), 
                        C = matrix(1, J, J),
                        tau = rep(0, J), scheme = "factorial")
                        

gcca.with.rgcca = rgcca(A= list(X,Y,Z), 
                        
                        C = C,
                        tau = rep(0, J), scheme = "factorial")



# X1 = Block1, ..., XJ = BlockJ, X_{J+1} = [X1, ..., XJ] # (J+1)*(J+1) Design matrix C
# C = matrix(c(0, 0, 0, ..., 0, 1,
#              0, 0, 0, ..., 0, 1,
#              0, 0, 0, ..., 0, 1,
#              ...
#              1, 1, 1, ..., 1, 0), J+1, J+1)
# Shrinkage parameters tau = c(tau1, ...,  tauJ, tau_{J+1})
# gcca.with.rgcca = rgcca(A= list(X1, ..., XJ, cbind(X1, ..., XJ)), C = C,
#                         tau = rep(0, J+1), scheme = "factorial")

#regCCA
data=list(X,Y,Z)
aaa=regCCA(data,reg=0)

# 2.5.10 SSQCOR
J=3
ssqcor.with.rgcca = rgcca(A= list(X,Y,Z), C = matrix(1, J, J),
                          tau = rep(0, J), scheme = "factorial")


J=3
ssqcor.with.rgcca = rgcca(A= list(X,Y,Z), C = matrix(1, J, J),
                          tau = rep(0, J), scheme = "centroid")

# Shrinkage parameters tau = c(tau1, ..., tauJ)
C = matrix(1, J, J) ; diag(C) = 0
maxbetb.with.rgcca = rgcca(A= list(X,W,Z),
C = C,
tau = rep(1, J), scheme = "factorial")



